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Abstract 


This  report  finishes  the  work  supported  by  the  Office  of  Naval  Research  as  cited  in 
the  acknowledgement.  In  particular  it  indicates  the  numerical  and  practical  feasibility 
of  the  tableau  approach  to  multifrequency  fault  diagnosis  previously  introduced  in 
reference  [l].  The  first  chapter  presents  a  26  component  fault  diagnosis  example 
based  on  a  linearized  model  of  a  video  amplifier  circuit.  The  diagnosis  of  the  example 
circuit  depends  on  the  solution  of  the  Tableau  Fault  Diagnosis  Equations.  The 
development  of  the  equations  as  well  as  a  solution  technique  are  the  subject  of  an 
earlier  report  (see  reference  [1]).  Chapter  2  investigates  the  properties  of  the  tableau 
fault  diagnosis  equations  under  the  assumption  that  the  number  of  simultaneous 
faults  in  a  system  under  test  is  limited.  The  approach  employs  multifrequency  testing 
to  extract  a  maximum  of  information  from  a  given  set  of  test  points.  The  analysis 
specifically  addresses  the  case  in  which  the  number  of  faulty  components  exceeds  the 
number  of  measurement  points.  The  development  includes  a  detailed  description  of 
the  solution  algorithm.  Several  example  problems  and  solutions  conclude  the  chapter. 
These  include  the  video  amplifier  circuit  which  is  used  to  illustrate  the  feasibility  of 
the  approach.  Appendices  are  included  which  detail  the  Fortran  code  of  the  computer 
programs  which  implement  the  fault  diagnosis  techniques  for  the  examples  of  the 
report. 


Chapter  1 

Completion  of  "Full  Diagnosis"  Research 


1.  Introduction 

The  purpose  of  this  report  is  the  completion  of  the  work  presented  in  [l].  This 
research  concentrated  on  the  "full  diagnosis"  problem,  i.e.  the  problem  of  determining 
system /circuit  parameters  from  multifrequency  output  measurements  under  the 
assumption  that  all  components  may  be  faulty  simultaneously.  Before  discussing  new 
material  we  summarize  the  major  contributions  of  the  previous  report  [  1  ]. 

First  is  the  development  of  a  set  of  Tableau  Fault  Diagnosis  Equations  based  on 
the  use  of  the  Component  Connection  Model  (CCM).  Due  to  the  nature  of  the  CCM  the 
fault  diagnosis  equations  and  the  associated  Jacobian,  utilized  in  the  solution 
algorithm,  have  an  elegant  and  sparse  structure.  Another  important  characteristic  of 
these  equations  is  that  in  many  cases  they  are  quadratic.  The  formulation  of  the 
diagnosis  equations  appears  in  summary  in  Chapter  2  of  this  report.  The  second 
contribution  of  [l]  is  the  presentation  of  a  theory  of  diagnosability  for  the  fault 
diagnosis  equations.  This  includes  a  test  for  diagnosability  based  on  the  computation 
of  the  rank  of  the  sparse  Jacobian.  Third  is  the  development  of  a  solution  algorithm 
for  the  nonlinear  Tableau  Fault  Diagnosis  Equations.  Of  particular  interest  is  a 
modification  which  exploits  the  quadratic  nature  of  the  equations  and  substantially 
improves  the  convergence  properties  of  the  solution  procedure.  Finally  [l]  presents 
several  illustrative  examples  which  include  algorithm  performance  data. 

The  third  section  of  this  chapter  completes  our  work  in  the  application  of  the 
fault  diagnosis  equations  to  the  "full  diagnosis"  problem  with  an  example  which 
illustrates  the  feasibility  of  the  Tableau  approach  for  large  systems.  The  second 
chapter  then  adapts  this  work  to  the  assumption  that  the  number  of  simultaneous 
faults  is  limited.  Before  proceeding  to  new  material  we  present  a  review  of  the  CCM 
and  the  derivation  of  the  fault  diagnosis  equations. 

2.  The  Component  Connection  Model  and  the  Fault  Diagnosis  Equations. 

Let  a  linear  system  have  N  components  where  the  i-th  component  is 
characterized  by  the  transfer  function  Zj(s,rj)  (s  is  the  Laplace  Transform  variable  and 
r*  is  a  parameter  which  characterizes  the  component).  Denote  the  i-th  component 
input  and  output  as  iij(s)  and  bj(s)  respectively.  Then 
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b,(s)  =  Z,(s.rj)a1(s)  (2.1) 

is  the  component  input/output  equation.  Next  define  the  composite  component 

input/output  vectors: 

a(s)  =  col  (ai(s),a8(s) . aN(s))  (2.2) 

b(s)  =  col  (bi(s),ba(s) . bN(s))  (2.3) 

and  the  composite  component  transfer  function 

Z(s,r)  =  block  diag(...1Zi(s1ri),...)  (2.4) 

where  r  =  col  (r!,r8 . rn).  The  composite  component  input  and  output  vectors  are 

related  by 

b(s)  =  Z(s,r)a(s)  (2.5) 

The  connection  laws  (  e.g.  KVL,  KCL)  are  then  expressed  in  terms  of  the  following 
equations: 

a(s)  =  Liib(s)+L18u(s)  (2.6) 


y(s)  =  L81b(s)+L22u(s) 


(2.7) 


where  u(s)  and  y(s)  are  the  circuit /system  input  and  output  vectors  respectively  and 
the  Ljj  are  determined  by  the  connections.  Equations  2.5,  2.6  and  2.7  form  the  CCM 
equations  [5]  and  have  a  frequency  domain  tableau  formulation: 


Z(s.r) 

-I 


0 

-Lizu(s) 


(2.8a) 


y(s)  =  L81b(s)+L2Zu(s) 


(2.8b) 


Suppose  the  circuit/system  modeled  by  equation  2.8  is  to  be  diagnosed:  that  is 
the  value  of  r  in  Z(sj,r)  is  sought.  Test  inputs  are  applied  at  q  different  test 
input/frequency  combinations  and  the  corresponding  outputs  are  measured.  Let  u(s() 
be  the  test  input  vector  and  yH(Sj)  be  the  test  output  vector,  for  i=l,2,...,q.  Because 
the  circuit/system  is  linear  all  components  of  u(sj)  and  yM(Sj)  are  phasors.  It  is 
possible  to  construct  a  set  of  fault  diagnosis  equations  of  the  form  [l]  : 


[Z(sj,r)  |  — V] 


LuVoi+a0(s,) 

Oj 


b0(si) 


(2.9) 


for  i  =  l,2 . q  ,  where 

(1)  q  is  thorn  of  '  input/frequency  combinations, 


(2)  b0(Si)  =  Lgfty^sO-Laaufsi)], 

(3)  La"/*  is  any  right  inverse  of  L2i, 

(4)  a0(Sj)  =  Lnb0(Si)  +  L^uCs,). 

(5)  V  is  a  matrix  whose  columns  span  the  null  space  of  L^, 

(6)  r  is  the  unknown  parameter  vector,  and 

(7)  Oj  is  a  vector  of  auxiliary  unknowns  which  characterizes  the  ambiguity  in  the 
solution  for  r  at  any  single  frequency,  sj. 

In  particular  if 

fi(r)  4  [Z(sj,r)  |  -V]  (2.10) 


gi(oj)  = 


LuVoi+a0(si) 
!  -£i 


(2.11) 


ft  ^  b0(Si) 

X  k  xcol  . otq.r] 

then  the  fault  diagnosis  equations  have  the  equivalent  form 

fi(r)gi(a,)-0i 

f(x)  =  ;  =© 

fq(r)s  q(£Lq)~0q 


(2.12) 

(2.13) 


(2.14) 


where  0  is  the  zero  vector.  Using  a  Newton-Raphson  scheme,  one  iteratively  solves[l] 
for  the  solution,  say  x*.  via 

JF(xk)[xk+1-xk]  =  -F(xk)  (2.15) 

where  xk  is  the  k-th  estimate  of  the  solution  to  equation  2.14  and  JF(.)  is  the  Jacobian 
of  F(.)  .  Because  of  the  product  structure  of  equation  2.13  and  the  sparsity  inherent  in 
2.9,  the  Jacobian  is  both  elegant  and  sparse.  Specifically 

dg,  af, 

0  •  aT-O-^Ca,) 

0rs  dfp 

W-^i2z)  ■  .  63(22) 

0.0 


JF(x)  = 


0 


(2.16) 
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Table  2. 

Nonzero  entries  of  L1Z  (26x4). 


row.column 

value 

1.1 

1 

9,2 

1 

19,3 

1 

22,4 

1 

Table  3. 

Nonzero  entries  of  La!  (6x28). 


row.column 

value 

row.column 

value 

1.1 

1 

2,2 

1 

2,4 

1 

2,6 

1 

3,9 

1 

4,19 

1 

4,26 

1 

5,22 

1 

6,5 

1 

6.6 

-1 

6,8 

-1 

6.10 

-1 

6,12 

-1 

6,17 

1 

6,18 

-1 

6,24 

1 

6,25 

-1 

The  computation  of  and  V  was  performed  via  the  IMSL  routines  LG1NF  and 
LSVDF  respectively  [8],  These  too  are  spare  matrices  with  their  nonzero  entries  given 
in  Tables  4  and  5. 
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Table  4. 

Nonzero  entries  of  L£iR  (26x6). 


row, column 

value 

row.column 

value 

1.1 

1 

2,2 

0.348154 

2,6 

0.0384615 

4,2 

0.346154 

4,6 

0.0384615 

5,2 

0.0364615 

5,6 

0.115385 

6,2 

0.307692 

6,6 

-0.0769231 

8,2 

-0.0384615 

8,6 

-0.115385 

9,3 

1 

10,2 

-0.0384615 

10,6 

-0.115385 

12,2 

-0.0384615 

12,6 

-0.115385 

17,2 

0.0384615 

17,8 

0.115385 

18,2 

-0.0384615 

18,6 

-0.115385 

19,4 

0.5 

22,5 

1 

24,2 

0.0384615 

24,6 

0.115385 

25,2 

26,4 

-0.0384815 

0.5 

25,8 

-0.115385 

Table  5. 

Nonzero  entries  of  V  (26x20) 


row, column 

value 

row, column 

value 

row, column 

value 

2,1 

0.808608 

3,11 

1 

4,1 

-0.428088 

4.2 

-0. 105202 

4,3 

-0.105202 

4,4 

-0.105202 

4,5 

0.0743889 

4,6 

-0.105202 

4.7 

0.105202 

4,8 

-0.105202 

4.9 

-0.62699 

4,10 

0.0743889 

5,1 

-0.0475652 

5,2 

0.109844 

5,3 

0.109844 

5,4 

0.109844 

5,5 

0.629435 

5,6 

0.109844 

5,7 

-0.109844 

5,8 

0.109844 

5,9 

0. 132594 

5,10 

0.629435 

6,1 

-0.380521 

6,2 

0. 105202 

8.3 

0. 105202 

6,4 

0.105202 

6,5 

-0.0743889 

8,6 

0.105202 

6,7 

-0.105202 

6.8 

0. 105202 

6.8 

0.62699 

8,10 

-0.0743889 

7,12 

1 

8.1 

0.0475652 

8,2 

0.890158 

8,3 

-0.109844 

8,4 

-0. 109844 

6,5 

0.0776715 

8.6 

-0.109844 

8.7 

0.109844 

8,8 

-0.100844 

8.9 

-0.132594 

8,10 

0.0776715 

10.1 

0.0475652 

10,2 

-0.109844 

10,3 

-0.109844 

10,4 

0.890156 

10,5 

0.0778715 

10,6 

-0. 109844 

10,7 

0.109644 

10,8 

-0. 109844 

10.9 

-0.132594 

10,10 

0.0770715 

11,13 

1 

12,1 

0.0475852 

12,2 

-0.109844 

12,3 

-0. 109844 

12,4 

-0.109844 

12,5 

0.0778715 

12,6 

0.890156 

12,7 

0.109844 

12,8 

-0.109844 

12,9 

-0.132594 

12,10 

0.0776715 

13,14 

1 

14,15 

1 

15,18 

1 

16,17 

1 

17,1 

-0.0475652 

17,2 

0.330293 

17,3 

0.336293 

17,4 

0.336293 

17,5 

-0.237795 

17,6 

0.336293 

17,7 

-0.336293 

17.8 

0.338293 

17,9 

-0.301105 

17,10 

-0.237795 

18,1 

0.0475652 

18,2 

-0.109844 

IB, 3 

0.890158 

18,4 

-0. 109844 

18,5 

0.0776715 

18,6 

-0.109844 

18,7 

0.109844 

18,8 

-0.109844 

18,9 

-0.132594 

18,10 

0.0778715 

19,5 

0.5 

19,10 

-0.5 

20,18 

1 

21,19 

1 

23,20 

1 

24,1 

-0.0475852 

24,2 

0.109B44 

24,3 

0.109844 

24.4 

0. 109844 

24,5 

-0.0776715 

24,6 

0. 109844 

24,7 

0.890156 

24,8 

0.109844 

24,9 

0. 132594 

24,10 

-0.0778715 

25,1 

0.0475652 

25,2 

-0.109844 

25,3 

-0.109844 

25,4 

-0.109844 

25,5 

0.0770715 

25,8 

-0. 109844 

25,7 

0.109844 

25,8 

0.890158 

25.9 

28.10 

-0.132594 

0.5 

25,10 

0.0770715 

26,5 

-0.5 

The  component  transfer  functions  are:  Z|(s,rj)  =  rjs  for  i=4, 5, 10, 11, 16, 17, 23, 24  and 
Zi(s.ri)  =  T|  for  the  remaining  i. 


The  nominal  component  values  appear  in  Table  8  which  includes  the  solution  results. 
Note  that  the  component  values  arc  scaled  to  improve  the  numerical  condition  of  the 
problem.  The  impedance  scale  factor  is  103  and  the  frequency  scale  factor  is  107. 
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For-  this  example  M=N=26  and  p=20.  Using  theorem  6.1  the  minimum  number  of 
input / frequency  combinations  is  q=3.  For  this  value  of  q  we  found  no  set  of  inputs  and 
frequencies  for  which  the  circuit  would  meet  the  diagnosability  condition  of  Theorem 

.2.  Increasing  q  to  4  we  And  that  the  following  inputs  and  frequencies  render  the 
circuit  diagnosable: 


u(sO  =  u(j.05)  = 


(7  la) 


u(s2)  =  u(j.2)  = 


(7.1b) 


u(s3)  =  u(j.l)  = 


(7.1c) 


u(s4)  =  u(jl.)  = 


(7.  Id) 


The  nominal  values  for  the  Oi  to  be  used  along  with  the  nominal  parameters  as  the  first 
solution  estimate  are  in  Tables  6  and  7. 
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Table  6. 

Nominal  values  for  oti  and  o^. 


a 

real  part 

Li 

imag.  part 

— .38813E— 01 

— .10578E— 01 

.10B73E+00 

.11 109E— 01 

.23505E+00 

— .23528E+00 

.25592E+00 

.20359E+00 

.90589E+00 

— .17871E+00 

.84608E  +00 

—.801 19E+00 

.31807E+00 

— .B8404E— 01 

.78337E  —01 

— .74524E— 01 

—.5 1649E-01 

.18481E+00 

— .10085E-01 

.84373E+00 

.29313E+01 

—  42463E+01 

.  12641E+01 

— .19982E+01 

.49954E+00 

.31802E+00 

.31785E+00 

-.23551  E+00 

.13377E— 01 

— .25264E— 01 

.39654 E+00 

— .44345E+00 

.63160E— 01 

.33443E— 01 

.44210E+00 

— .77873E— 01 

— .22943E— 02 

— .91910E— 02 

.22977E— 01 

— .57358E— 02 

.35685E-01 
.79972E— 01 
.  15327E+00 
.  18606E+00 
.37B58E+00 
.26385E+00 

•  23297E— 01 
.10206E+00 

•  54178E— 01 
.1507BE+00 
.B3646E+00 
.35621 E+00 
.  17659E+00 
.  16082E+00 
.69624E— 02 
.33091  E+00 
■  94932E— 01 
.13150E+00 
.  1B614E— 02 
.24575E-01 


— .22477E— 02 
— .  16134E— 01 
— .11376E+00 
— .36613E+00 
-.28573E— 01 
.2S066E+00 
— .28580E— 01 
— .43404E— 01 
.39788E— 01 
•  52261E— 01 
— .55977E— 01 
-.  17659E+00 
.35621  E+00 
-.29919E-01 
-  94932E— 02 
— .521 13E— 01 
.69824E-01 
— .29041E— 01 
— .24575E— 02 
.  1B614E— 01 
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Table  7. 

Nominal  values  for  otg  and  a*. 


a 

real  part 

8 

imag.  part 

a 

real  part 

4 

imag.  part 

•21705E— 01 

.83072E— 03 

— .37897E-01 

— .54919E— 03 

— .14510E+00 

-.54401 E-01 

.71877E— 01 

— .54767E— 02 

.2591  BE +00 

— .21989E+00 

.77880E— 01 

— .27958E— 01 

— .18065E  +00 

— .18939E+00 

•  17227E+00 

.88386E— 01 

•  12023E+01 

—  .33955E+00 

—  .52193E+00 

— .52442E— 03 

— .37128E+00 

■20523E— 01 

.90419E— 01 

-.25391  E-01 

.20897E  +00 

.82907E— 01 

.B0488E+00 

.28640E+00 

— .17104E+00 

— .54508E-01 

.31594E— 01 

•  17439E+00 

.1 1877E+00 

.10451E+00 

.  10015E+00 

.72395E— 02 

— .99168E+00 

.42984E+00 

■  34497E+00 

-.905375-03 

— .77921E— 01 

— .93854E— 01 

.53431E+00 

•  44778E— 01 

.251 16E+00 

— .79270E— 01 

— .  19795E— 01 

■30812E— 01 

.39635E-01 

.12558E+00 

— .  1 54 06 E +00 

— .98974E— 01 

-.32688E+00 

— .42537E— 01 

.88285E— 01 

.72916E— 02 

.40200E  —01 

— .16473E— 01 

■  59702E— 03 

-.22681 E-02 

.97877E+00 

.89241  E-01 

.23887E+00 

.  13B87E— 01 

.82365E— 01 

.20100E+00 

.11340E+00 

.29851  E-01 

.85552E— 01 

.13844E— 01 

— .22188E— 01 

.02441  E-01 

— .28225E— 02 

.84917E— 04 

— .40096E— 02 

•  17967E— 01 

— .32459E— 03 

—.141 12E— 01 

— .89B34E+00 

— .2004BE+00 

We  next  establish  a  set  of  "actual"  parameter  values  which  are  to  be  determined 
from  the  measurement  data.  These  values  likewise  appear  in  Table  B.  The  simulated 
measurements  corresponding  to  the  actual  parameters  and  the  inputs  of  equation  7.1 
are: 


y“(j05) 


■  30131E+01 

— .41295E+01 

.54763E— 01 

— .74921E— 01 

.13175E  +  Q1 

+j 

— .10620E+01 

.  10894E+01 

35500E+00 

.23035E-01 

— .77400E  +  00 

.36800E+00 

.7B442E  +  00 

(7.2a) 


yM(j.2)  = 


yM(j.l)  = 


.76267E+00 
.  13B62E-01 
.79305E  +  00 
.37143E+00 
.42976E+00 
.66459E+00, 


.  10196E+00 
-.1 050-1 E -02 
— .  10B65E+00 
.  17885E+01 
-.26883E+00 
,— .45267E  +  00, 


.30754E— 01 
—  .54854E— 03 
— .  14179E+00 
-.64531E-02 
-.40539E-01 
.41443E— 01, 


—  .94506E— 01 
-.17106E-02 
— .78B25E— 01 

—  .75712E  +  00 
— .  1 1 197E+00 
.-.61275E-02, 


(7.2b) 


(7.2c) 
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yM(ji) 


.50681E+00 

.38099E— 01 

9210BE-02 

.69435E— 03 

.50637E+00 

+i 

.27279E— 01 

.22685E+00 

.35169E-02 

54900E+00 

— .17301E+00 

,74654E+00, 

•  28013E-02, 

(7.2  d) 


Table  B  summarizes  the  results  of  the  solution  algorithms.  The  first  solution  used 
the  usual  Newton-Raphson  iteration  step  while  the  second  solution  used  the  modified 
algorithm  discussed  in  Chapter  5  of  [1].  Both  solutions  employed  a  FORTRAN  program 
compiled  and  executed  on  a  VAX-1 1/7B0  with  single  precision  arithmetic.  Use  of  the 
usual  or  modified  iteration  is  a  program  option.  This  program  is  included  as  Appendix 
A  of  this  report.  Although  both  solution  methods  converged  for  this  example  the 
modified  algorithm  required  7  iterations  (17  min.)  to  converge  compared  to  17 
iterations  (36.8  min.)  for  the  Newton-Raphson  algorithm.  The  execute  times  for  the 
programs  are  somewhat  long  due  to  the  fact  that  sparse  matrix  techniques  were  not 
used.  The  Jacobian  for  this  example  is  a  208x186  matrix  with  5.4%  of  its  elements 
nonzero. 
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a 


s 

w  ■ 


k. 


\r 


Table  8. 

Parameter  values  and  solution  results. 


component 

units 

nominal 

actual 

solution  1 

solution  2 

1 

ohm 

12 

11 

n 

11 

2 

mho 

0.1 

0.11 

0.109 

0.109 

3 

ohm 

56.7 

55 

55 

55 

4 

farad 

50 

100 

100.2 

100.2 

5 

farad 

5 

4.5 

4.5 

4.5 

6 

mho 

10 

9 

9.016 

9.017 

7 

ohm 

30 

37 

37 

37 

a 

mho 

0.1 

0.15 

0.1498 

0.1498 

e 

ohm 

10 

a 

a 

8 

10 

farad 

50 

65 

05 

65 

11 

farad 

5 

6.2 

6.2 

6.2 

12 

mho 

10 

14 

14 

14 

13 

mho 

0.3 

0.22 

0.22 

0.22 

14 

ohm 

10 

11 

11.03 

11.03 

15 

ohm 

2 

2.5 

2.5 

2.5 

16 

farad 

60 

45 

45.04 

45.04 

17 

farad 

5 

4.6 

4.8 

4.8 

ia 

mho 

10 

9 

9.007 

9.007 

19 

ohm 

10 

11 

11 

11 

20 

mho 

0.3 

0.33 

0.33 

0.33 

21 

ohm 

10 

11 

11.23 

11.23 

22 

ohm 

10 

12 

12 

12 

23 

farad 

50 

49 

49 

49 

24 

farad 

5 

6 

6 

6 

25 

mho 

10 

3.8 

3.798 

3.798 

26 

ohm 

0.78 

0.7 

0.7 

0.7 

4.  Summary 

The  preceding  example  completes  the  Tableau  Fault  Diagnosis  Equations 
documented  in  [1]  by  demonstrating  the  feasibility  of  their  use  for  a  large  system. 
These  equations  are  easily  computed,  requiring  no  matrix  inversions  in  the  their 
construction  or  in  the  construction  of  the  associated  Jacobian.  The  polynomial  order 
of  the  Tableau  Fault  Diagnosis  Equations  is  dependent  on  the  component 
characteristics  and  not  the  size  of  the  system.  For  many  systems  this  means  that  the 
equations  are  quadratic,  a  fact  which  has  been  exploited  to  improve  the  convergence 
properties  of  the  solution  algorithm.  Finally  a  substantial  improvement  in  the 
efficiency  of  the  solution  algorithm  is  possible  upon  the  application  of  sparse  matrix 
techniques. 


Chapter  2 

Fault  Diagnosis  in  the  Tableau  Context 
with  the  Assumption  of  Limited  Simultaneous  Faults 


1.  Introduction 

The  discussion  of  the  previous  report  [1]  presumed  that  ALL  parameters  of  a 
circuit/system  could  differ  significantly  from  their  nominal  values.  This  assumption  is 
appropriate  when  there  is  some  form  of  interdependence  among  the  parameters.  An 
example  of  this  is  the  parameter  set  of  a  transistor  model.  A  transistor  functioning 
abnormally  could  cause  all  the  parameters  of  a  linearized  model  to  change 
significantly.  However  parameter  failures  in  many  circuits/systems  are  often 
statistically  independent.  In  this  case  the  likelihood  of  more  than  a  few  simultaneous 
failures  is  extremely  small.  The  fact  that  most  parameters  are  at  or  near  their 
nominal  values  represents  information  useful  to  the  solution  process  and  leads  to  the 
following  objective  for  this  chapter: 

Utilize  the  Tableau  fault  diagnosis  equations  developed  in 
[l]  to  solve  for  the  circuit/system  parameter  vector,  r, 
given  i)  the  measurement  data  and  ii)  the  assumption  that 
at  most  nf  of  the  N  parameters  differ  from  nominal  where 
the  integer,  nf,  denotes  the  "number  of  assumed  faults". 

Several  compelling  motivations  underly  this  recasting  of  the  fault  diagnosis 
problem.  First  it  simplifies  the  calculations  since  the  solution  algorithm  deals  with 
fewer  equations  and  unknowns.  In  fact  there  are  recent  well  documented 
circumstances  [2-4]  in  which  the  computations  proceed  via  linear  methods.  More 
important  than  the  simplification  of  the  diagnosis  equations  is  the  fact  that  the  fault 
diagnosis  process  requires  fewer  test  points  when  the  assumption  of  a  limited  number 
of  faults  is  valid.  Call  the  number  of  test  points  (outputs)  rv  Each  of  the  n«,  test 
points  is  a  source  of  information  which  in  composite  should  suffice  to  determine  the 
parameters.  Intuitively,  the  amount  of  information  necessary  is  proportional  to  the 
dimension  of  the  parameter  space.  Since  the  limited  fault  assumption  forces  the 
solution  to  lie  in  a  lower  dimension  subspace  of  the  parameter  space,  it  follows  that  its 
use  in  the  diagnosis  process  should  require  fewer  test  points  than  a  complete 
diagnosis  (i.c.  no  limited  parameter  failure  assumption).  Keep  in  mind  however  that  a 
tradeoff  arises  between  the  reduction  of  test  points  and  the  ease  of  solving  the 
resulting  equations.  For  a  given  value  of  nf  as  the  number  of  test  points,  n0,  decreases 
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the  difficulty  of  solving  the  fault  diagnosis  equations  should  increase. 

Current  research  in  the  literature  has  concentrated  on  simplifying  the  fault 
diagnosis  calculations,  For  example  Huang,  Lin  and  Liu  [2]  have  developed  a  node¬ 
fault  approach  which  solves  for  the  change  in  the  branch  admittances  of  a  linear 
network  using  voltage  measurements  at  a  set  of  "accessible”  nodes  in  the  network.  If 
the  number  of  accessible  nodes  is  m+1  (including  the  reference  node)  and  the  number 
of  faults  is  k,  this  approach  requires  that  k<m  for  the  fault  to  be  uniquely  determined 
using  linear  techniques.  (Note:  For  this  approach  the  number  m  is  analogous  to  the 
number  n0  in  the  Tableau  approach  and  k  is  analogous  to  nf.)  For  any  set  of  k  faulty 
branches  where  k<m,  the  set  of  possible  values  for  the  network  response  lies  in  a  k- 
dimensional  linear  subspace  of  Rm.  If  a  measurement  does  not  lie  in  any  such 
subspace  then  k^m  and  it  is  not  possible  to  determine  the  fault  uniquely.  If  on  the 
other  band  the  measurement  does  lie  in  a  subspace  correspdnding  to  a  particular  k- 
fault  then  the  branch  admittances  corresponding  to  that  fault  are  solvable. 

Likewise,  Saeks  et  al.  [3]  have  introduced  a  fault  diagnosis  algorithm  in  the 
limited  failure  context  based  on  the  CCM.  The  authors  recognized  the  excessive 
computational  cost  associated  with  solving  a  set  of  fault  diagnosis  equations  for  every 
possible  fault  combination  and  their  approach  circumvents  this  problem.  Their 
algorithm,  applicable  to  nonlinear  as  well  as  linear  networks,  has  the  following 
structure: 

i)  Partition  the  components  into  two  groups.  Assume  that  the 
group  1  components  are  "good".  Using  the  measurement 
data  and  the  characteristics  of  the  group  1  components 
determine  the  inputs  and  outputs  of  the  group  2 
components  which  would  give  rise  to  the  measurement 
data. 

ii)  Test  each  component  in  group  2  by  determining  if  the 
component  outputs  and  inputs  are  consistent  with  the 
nominal  component  characteristic.  If  the  inputs  and 
outputs  of  component  i  are  consistent  with  its  nominal 
characteristic  then  component  i  is  assumed  good,  if  not 
then  no  decision  is  possible. 

iii)  Repartition  the  components  including  all  components 
found  to  be  good  in  group  1.  Go  to  step  ii  and  continue  the 
process  until  group  1  consists  entirely  of  good  components. 

Several  remarks  concerning  this  procedure  are  relevant  to  the  development  of 
the  algorithm  in  this  paper.  First  the  decision  rule  used  in  [3]  to  establish  the  identity 
of  the  good  components  is  exact  for  single  component  failures  but  for  multiple  faults 
the  decision  rule  must  be  based  on  the  assumption  that  the  effects  of  multiple  faults 
cannot  cancel.  This  assumption  is  reasonable  when  the  "good"  components  are 
represented  by  parameters  which  are  exactly  nominal  but  may  break  down  when 
actual  values  of  the  good  parameters  arc  randomly  spread  around  nominal  due  to 
production  variations.  Second,  the  algorithm  in  [3)  requires  that  the  number  of  Lest 
points  for  a  linear  system  be  sufficient  to  permit  the  compulation  of  the  inputs  and 
outputs  of  the  components  under  test  with  single  frequency  lest  data.  Third,  use  of 
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the  multiple  fault  decision  rule  requires  that  the  number  of  faults  be  strictly  less  them 
the  number  of  components  under  test.  (Note:  These  two  requirements  are  equivalent 
to  the  restriction  in  the  Huang,  Lin  and  Liu  approach  [2]  that  k<m.)  Finally  although 
this  method  employs  the  CCM  it  differs  fundamentally  from  the  Tableau  Approach 
presented  in  this  paper. 

A  Anal  example  of  the  use  of  the  limited  fault  assumption  is  that  of  Biemacki  and 
Bandler  who  developed  an  approach  to  multiple  fault  location  for  linear  networks. 
Here  the  faults  are  modeled  as  loads  applied  to  an  invariant  network  which  represents 
the  system  in  its  nominal  state  [4].  Their  fault  diagnosis  equations  build  on  voltage 
measurements  taken  at  a  single  test  frequency.  The  identification  of  the  faults 
depends  upon  checking  the  consistency  of  a  set  of  linear  equations.  Like  the  other 
approaches  the  linearity  of  the  fault  diagnosis  equations  requires  that  the  number  of 
voltages  measured  be  greater  than  the  number  of  simultaneous  faults  to  be  located. 
(Note:  As  for  the  previous  approach  this  assumption  is  equivalent  to  the  restriction  in 
the  Huang.  Lin  and  Liu  approach  [2]  that  k<m  or  equivalently  in  the  Tableau  context 
that  nf<n0.) 

All  the  above  approaches  share  the  constraint  that  nf<rio.  The  algorithm  for 
limited  fault  analysis  developed  in  this  chapter  does  not  require  that  the  number  of 
faults,  nf,  be  less  than  the  number  of  measurement  points,  tv  Philosophically,  the 
idea  is  to  use  the  multifrequency  techniques  developed  thus  far  to  "squeeze"  as  much 
information  as  possible  out  of  a  given  set  of  test  points.  The  next  section  presents 
several  simple  examples  in  the  context  of  the  Tableau  Fault  Diagnosis  Equations 
applied  to  the  limited  fault  assumption.  These  examples  will  clearly  illustrate  how  the 
relationship  between  and  n,  affects  the  approach  to  the  solution. 


2.  Motivational  Examples 


The  purpose  of  this  section  is  to  illustrate  the  properties  of  the  Tableau  Fault 
Diagnosis  Equations  under  the  assumption  that  the  number  of  parameters  which 
deviate  significantly  from  nominal  is  limited.  The  examples  are  designed  to  highlight 
the  differences  which  occur  as  the  relation  between  nf  and  no  changes  and  to  serve  as 
background  to  the  development  of  the  solution  algorithm. 

Consider  the  example  shown  in  Figure  2.  The  CCM  equations  for  this  example  are: 

bi(s) 
bz(s) 
b3(s) 
b4(s) 


ra/  s  0 

ai(s) 

a2(s) 

0  r3 

a3(s) 

r4 

a4(s) 

a>(s) 

az(s) 

a3(s) 

a4(s) 

-1 

0 

1 

0 


0 

-1 

0 

0 


bi(s) 

i 

b2(s) 

0 

b3(s) 

4* 

0 

b4(s) 

1. 

u(s) 


(2  lb) 


•  •  -w  / 


bi(s) 


r  nl  fbi(s) 

10  0  1  bg(s)  fo|  ,  v 
0  10  0]  b3(s)  +  [oj  u(s' 

Ms) 


(2.1c) 


where 


r2 

r  ~ 

r3 


Gi 

1/C2 

G3 

g4 


i(s)  _  '  I 

2(s)  "  k 


The  right  inverse  and  null  space  basis  for  L^i  are: 


1*21  “ 


1 

0 

0 

0 

and  V  = 

0 

1 

-1 

0 

The  nominal  value  for  the  parameter  vector,  r0  is  r0  =  col[l  111]  and  0^=2. 

Case  1:  nf=  1  (Component  1  faulty) 

.  Assume  that  only  one  parameter  value  differs  from  nominal  and  that 
measurements  occur  for  a  single  real  test  frequency,  s=Sj.  Recall  that  the  fault 
diagnosis  equation  is: 


[z(s!,r)  |  — V 


LuVot^aoCs!) 

0.1 


-b0(sj)  =  0 


where: 


i.  b0(si)  =  L1^[yH(s1)-L2Zu(si)]: 

ii.  a0(s!)  =  Li,b0(s,)+L12u(s1); 

iii.  The  columns  of  V  span  nullfLgi]; 

iv.  L£jR  is  any  right  inverse  of  L2]; 

v-  yM(si)  is  the  measurement  (an  n0  dimensional  vector): 

vi.  The  dimension  of  null[L2j]  is  p; 

vii.  a,  =  col(a1(s,),a2(si),...,ap(s1)). 

Next,  rearrange  the  fault  diagnosis  equation  2.2  to  produce: 


[v-ZCsj.rjLuV  a,  =  Z(s1,r)a0(si)-b0(s1) 


Substituting  the  information  for  this  example  into  equation  (2.3)  yields: 

1  °  1  f  ,  si  -1  “r> 

-r2/si  r2/s,  |r,Uo  1  I  r2/s,  -1 


iu(s,) 

r\ 


\rz/si 


•i  I  rv>.Vc^l 
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Notice  the  following  characteristics  of  this  set  of  equations: 

i.  There  are  4  equations  and  6  unknowns  (4  parameters  and  2  ambiguity 
variables); 

ii.  Each  rs  appears  only  in  the  i'th  equation; 

iii.  If  the  nominal  numerical  values  for  r  replace  the  variables,  the  equation 
set  is  linear  with  the  and  a2  as  unknowns; 

iv.  If  one  parameter  is  faulty  then  a  solution  for  at  and  a2  must  exist  upon 
deleting  the  equation  containing  the  nominal  value  of  the  faulty 
component. 

The  above  characteristics  suggest  a  solution  method.  Namely  eliminate  the  i’th 
equation.  If  the  remaining  equations  are  not  consistent  then  rj  is  not  faulty  (or  the 
single  fault  assumption  is  not  valid). 

To  illustrate  such  a  solution  method  let  component  1  be  faulty  with  rj  =  2.  Given 
the  test  frequency  st  =  1  and  input  u(Si)  =  1  the  corresponding  test  outputs  are: 

-  U] 


Equation  2.3  becomes: 

1 

-1 

0 

-1 


-1.5 

al 

1.5 

"3 

.5 

1.  . 

(2.6) 


y“(si) 

y“(si) 


where  =  a,(l)  and  a2  =  a2(l).  The  results  of  the  successive  elimination  of  each  or 
the  rows  of  equation  2.6  appear  in  summary  form  in  Table  9. 
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Table  9. 

Summary  of  Consistency  Check  (Component  1  Faulty). 


Row  Elim. 

Resulting  Equation 

Solution 

1 

1 

I 

1 

1 

2 

I 

1  0 

0  1 

-1  0 

1 

1 

1 

Inconsistent 

3 

1 

I 

1 

1 

1 

Inconsistent 

4 

1 

1 

1 

H 

1 

Inconsistent 

The  data  in  Table  9  clearly  indicates  that  the  only  possible  single  fault  for  the 
given  measurement  data  is  component  1.  To  determine  the  parameter  values  simply 
substitute  the  otj  values  into  2.7a  and  2.7b  to  compute  actual  component  inputs  and 
outputs: 

b(t)  *  b0(l)  +  V*x  (2.7a) 


a(l)  =  a0(l)  +  Lj  jVot , 


(2.7b) 


Then  use  the  composite  component  transfer  function  matrix  equation  b(l)  =  Z(l,r)a(l) 
to  compute  r.  For  this  example  equation  2.7  yields: 


b(l)  = 


'1 ' 

.5 

.5 

.5 

and  a(  1)  = 

.5 

.5 

,1 

,1 

and  the  parameter  values  are  ri  =  2  and  rz  =  r3  =  r4  =  1  as  expected. 


Case  2:  nf=  1  (Component  2  faulty) 


The  solution  to  a  limited  fault  problem  is  not  necessarily  unique.  To  demonstrate 
this  repeat  the  example  with  component  2  faulty.  Let  rz  =  2  with  the  remainder 
nominal.  The  test  outputs  for  this  case  given  s x  =  1  and  u(si)  =  1  are: 


yf*(si) 

1.6' 

yaH(s.) 

1.4  , 

(2.8) 


Substitute  into  equation  2.4  to  get: 
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1  0 

-1  1 

ai 

-1. 

1.2 

0  1 

a2 

.4 

-1  0 

1. 

Table  10  contains  the  results. 


(2.9) 


Table  10. 

Summary  of  Consistency  Check  (Component  2  Faulty). 


Row  Elim. 

Resulting  Equation 

Solution 

1 

-1  1 

0  1 

-1  0 

<*2 

= 

1.2 

.4 

1.  , 

Inconsistent 

2 

1  0 

0  1 

-1  0 

«z 

= 

-1.' 

.4 

.  1.  , 

az 

■  ul 

3 

1  0 

-1  1 

-1  0 

ai 

a2 

~ 

-1. 

1.2 

1. 

az 

4 

1  0 

-1  1 

0  1 

= 

-1.5 

1.5 

.5 

Inconsistent 

The  data  in  Table  10  indicates  that  there  are  two  possible  single  faults  which 
satisfy  the  fault  diagnosis  equations.  The  values  of  a  from  the  second  row  in  the  table 
correspond  to  parameter  values:  rz  =  2  with  the  remainder  nominal.  Those  of  the  third 
row  correspond  to  r3  =  .5  with  the  remainder  nominal.  This  illustrates  the  potential 
for  ambiguity  in  the  solutions.  It  is  easy  to  see  for  this  simple  circuit  that  the 
ambiguity  is  generic;  that  is,  it  will  always  occur  for  this  test  arrangement  whenever 
parameters  two  or  three  are  faulty.  Although  ambiguities  are  always  possible,  it  would 
be  convenient  to  be  able  to  avoid  such  generic  ambiguities.  A  later  section  will 
develop  a  test  to  determine  if  such  a  problem  exists  for  a  given  diagnosis  situation. 

To  summarize  cases  1  and  2  of  the  example  consider  the  following  characteristics. 
First,  the  number  of  faults,  nr,  is  strictly  less  than  the  number  of  test  points,  n„.  Next, 
each  fault  generated  a  new  set  of  equations  from  the  original  fault  diagnosis  equations. 
When  a  solution  failed  to  exists  for  the  set  corresponding  to  some  fault  then  that  fault 
was  not  possible.  The  existence  of  solutions  for  more  than  one  fault  possibility  may  be 
generic  to  the  circuit/system  test  point  combination.  Finally,  notice  that  the  solution 
process  proceeded  via  linear  methods.  This  is  possible  because  nj  <  tv  Here  the  use 
of  measurements  at  a  single  test  frequency  provides  sufficient  data  for  the  solution 
process  >n  that  the  set  of  measurements  for  each  fault  combination  must  lie  in  a 
linear  n.crpacc  of  the  entire  measurement  space.  Although  the  solution  to  these 
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cases  is  computationally  different  from  the  approaches  discussed  earlier  [2-4]  the  use 
of  linear  methods  and  its  underlying  justification  is  an  essential  commonality  which 
makes  them  all  philosophically  equivalent. 


Case  3:  nf=2  (Components  1  and  2  faulty) 

Consider  the  possibility  that  the  circuit  modeled  by  equation  2.1  has  two 
simultaneous  faults.  There  are  six  possible  fault  combinations,  i.e.:  rj,^  faulty;  ri.rs 
faulty;  etc.;  all  combinations  of  four  parameters  taken  two  at  a  time.  Consider 
equation  2.4  under  the  circumstance  of  two  simultaneous  faults.  As  before  there  are  4 
equations  and  6  unknowns,  Since  a  fault  combination  includes  two  parameters 
eliminating  the  equations  which  contain  an  assumed  fault  while  setting  the  remaining 
parameters  to  nominal  leaves  two  equations  in  two  unknowns.  This  means  that  a 
solution  for  ot]  and  a2  is  likely  to  exist  for  each  possible  fault  combination. 


Suppose  that  the  actual  fault  combination  for  this  example  is  components  1  and  2 
(r!  =  2  and  ra  =  2)  and  that  st  =  1.  The  test  outputs  for  this  case  are: 


y“(si) 

_  ri.86 

yl'(si) 

“  1.57 

(2.10) 


As  before  use  the  nominal  parameter  values  in  equation  2.4  with  the  understanding 
that  two  of  the  four  equations  must  be  incorrect  since  two  parameters  are  faulty. 
Equation  2.4  becomes: 


-1 

-1 

l' 

-1.43 

«i 

0 

1 

1 

-1 

[1.88]  _ 

1.29 

az 

0 

T 

0 

1 

1.57  J  ' 

.57 

,1 

0 

0 

,  1. 

(2-11) 


Table  11  summarizes  the  computations  of  the  a(  and  the  corresponding  r  for  each 
possible  fault. 


Table  11. 

Summary  of  Parameter  Computation  at  Sj  =  1  (Components  1  and  2  Faulty). 


Rows  Elim. 

Resulting  Equation 

Solution 

Parameters 

1&2 

0  1 
— !  0 

«i  _  r.5?i 
<*z  L 1.  J 

— 

<*! 

l  =  M 

1.57J 

r,=2  r3=  1 
rz=2  r4=l 

1&3 

1 

■ 

«i  _  fl.29 

a2  ~  l  1-  . 

1 

1 

1 

-  f-l1 

“  L29j 

r  i=2  r3=,5 
rz=l  r4=l 

1&4 

1 

M 

811 

1 

1 

1 

n 

■ 

rt  =  2.7  r3=l 
r2=  1  r4=.72 

2&3 

1 

■ 

1 

Inconsistent 

Inconsistent 

1 

99 

1 

I 

_  f— 1.43 
-  1  .57 

■ 

r,=l  r3=l 
r2=  -4  r4=  1 .43 

3&4 

n 

■ 

91 

I 

1 

1 

m 

■ 

r,=  1  r3=  -.25 
rz=  1  r4=  1.43 

As  expected  the  computation  produces  inconclusive  results.  Clearly  the 
resolution  of  the  ambiguity  in  the  data  of  Table  11  requires  additional  information.  To 
acquire  the  needed  information  repeat  the  same  procedure  using  a  second  test 
frequency  (the  verification  step).  Let  sz  =  2  for  which  measurement  at  the  output 
would  yield: 


y“(s2) 

_  fe. 

yi*(sz) 

~  l5j 

(2.12) 


As  before  substitute  this  information  into  equation  2.4  to  obtain: 


ai 

1 

0 

-1  -1 

.5  -1 

f2. 

-1.5 

.5 

0 

+ 

0 

1 

1.5 

.5 

,1 

0 

0 

,  1. 

(2.13) 


Table  12  contains  the  summary  of  results  for  the  computation  at  this  test  frequency. 
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Table  12. 

Summary  of  Parameter  Computation  at  st  =  2  (Components  1  and  2  Faulty). 


Rows  Elim. 

Resulting  Equation 

Solution 

Parameters 

1&2 

■ 

m 

1 

I 

■ 

r,=2  r3=l 
rz=2  r4=l 

1&3 

1 

t 

IS 

■ 

1 

AS 

■ 

r,=2  ra=0 
ra=l  r4=l 

1&4 

-.5  .5 

0  1 

i 

1 

i 

ii 

"~T 

01  b. 

rt=3.  r3=  1 
rz=l  r4=.5 

2&3 

Si 

i  _  f— 1.5 
a  t  1 

Inconsistent 

Inconsistent 

2&4 

1 

_  r-i.5i 
“  1  .5  j 

1 

1 

■ 

ri=l  r3=  1 
r2=«  r4=1.5 

3&4 

H 

SI 

B 

■ 

1 

1 

3B1 

■ 

r,  =  l  r3=  — 1 
r2=l  r4=  1.5 

The  final  step  in  this  two-fault  case  is  the  comparison  of  the  Tables  11  and  12.  The 
only  common  solution  between  the  two  tables  is  that  for  the  component  1  and  2  fault 
combination.  This  fault  set  therefore  is  the  only  one  which  simultaneously  satisfies  the 
equations  derived  from  measurements  at  both  test  frequencies. 

This  case  for  which  the  number  of  faults,  nf,  equals  the  number  of  test  points,  n«j, 
prompts  the  following  observations:  In  general  a  single  real  test  frequency  is 
insufficient  to  resolve  the  ambiguity.  Data  from  the  measurement  at  a  second  test 
frequency  is  necessary.  Since  the  equations  resulting  from  the  different  test 
frequencies  are  solved  separately,  linear  techniques  are  still  possible.  The  only 
significant  difference  between  this  case  and  the  previous  two  lies  in  the  necessity  for 
the  verification  step  (the  second  test  frequency).  This  step  was  not  necessary  in  cases 
1  and  2  since  the  equations  used  to  solve  for  the  a’s  were  overspecified  due  to  the  fact 
that  nf  <  n0. 


Case  4:  nf=2,  n^l 

Finally  consider  the  same  example  (Figure  2)  with  one  test  point  instead  of  two. 
The  CCM  equations  are  the  same  as  equation  2. 1  except  that  the  output  equation  is: 


y,(s)  =  [1001] 


b,(s) 

ba(s) 

b3(s) 

b4(s) 


+  [0]  u(s) 


(2.14) 


Accordingly  the  right  inverse  and  null  space  basis  of  L21  change  to: 
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10  0 
0  1  0 
0  0  1 
-10  0 

The  fault  diagnosis  equations  at  a  single  test  frequency  St  for  u(si)  =  1  are: 

l-a2(Si) 

r *  0  0  0  -1  0  0  y“(si)  +  a1(s1)-ag(si)  u,  * 

0  rj/s,  0  0  0  -1  0  a2(s,)  ^ 

0  0  r3  0  0  0  -1  aiJSj)  “  0  =0  (215) 

0  0  0  r4  1  0  0  oi(s!)  0 

03(Si) 

For  the  previous  cases  no  &rif.  To  illustrate  the  properties  of  the  fault  diagnosis 
equations  for  nr  >  n<,  let  the  number  of  allowable  simultaneous  faults  be  two  (nj  =  2) 
and  notice  that  n,j  =  1.  Under  these  circumstances  observe  that  for  equation  2.15,  the 
elimination  of  two  equations  corresponding  to  an  assumed  fault  combination  with  the 
remaining  two  parameters  set  to  nominal  leaves  two  equations  with  three  unknowns. 
This  means  that  the  ambiguity  must  be  resolved  by  the  use  of  several  test  frequencies. 
The  resulting  equations  must  be  solved  simultaneously  with  the  unfortunate 
consequence  that  the  solution  method  cannot  employ  linear  techniques. 

This  case  represents  a  more  general  limited  fault  problem  than  the  previous 
cases  since  the  number  of  test  outputs,  n^,  may  be  smaller  than  the  number  of 
assumed  faults.  nf.  Since  this  more  general  approach  is  the  subject  of  the  remainder 
of  this  chapter  the  actual  solution  process  will  appear  in  a  later  section. 

Current  research  [2-4]  has  concentrated  on  the  fault  diagnosis  problem  with  the 
restriction  that  n0<nf.  A  motivation  for  this  is  that  the  solution  method  may  use  linear 
techniques  as  illustrated  by  all  but  the  last  case.  Unfortunately  the  restriction  to 
linear  methods  fails  to  exploit  all  of  the  information  available  at  the  test  points.  Thus 
problems  such  as  case  4  are  not  solvable  via  such  methods.  The  remainder  of  this 
chapter  considers  the  development  of  a  more  general  algorithm  for  limited  fault 
analysis  in  the  CCM  context.  Here  the  restriction,  n0<nf,  will  not  apply  The  method 
will  exploit  all  the  information  available  at  a  given  set  of  test  points  by  utilizing 
multifrequency  testing. 

3.  Introduction  to  Notation 

The  purpose  of  this  brief  section  is  to  develop  the  notation  to  assist  in  describing 
the  nf-fault  solution  algorithm  presented  in  the  next  section.  To  this  purpose  we 
define  the  following: 

<i,,i2 . inf>  £  The  fault  index 

It  is  an  ordered  nr— tuple  of  positive  integers  subject  to  the  following: 

>'i (n.i) 
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n(N,nf)  s  j  The  set  of  all  possible  rif-tuples  satisfying  3.1  { 

There  are  Cn:i1i  elements  in  fi(N,nf).  (Note:  CN:„f  denotes  the  number  of  combinations  of 
N  things  taken  nf  at  a  time.) 

If  y  =  <ii,i2 . inj>  then  let  M7  be  the  following  matrix  with  nf  columns: 

Mr*  [ei,|el8l  -|einf]  (3.2) 


where  e)  is  the  N-dimensional  unit  vector  with  a  1  in  the  i-th  position  and  0’s  elsewhere. 
For  the  fault  -yen(N,nf)  define  the  fault  space  P7  in  the  following  way: 

P7=  |r  |r  =  r0  +  My5,  peR^j  (3.3) 

We  say  that  the  fault  index  y  "represents"  the  parameter  r  if  rePr  In  other  words  y  is 
a  specific  nf-fault  combination. 

Note:  Limiting  the  number  of  allowable  faults  to  nf  means  that  there  are  at  most  nf 
faults  since  any  fault  space  with  less  than  nf  faults  is  contained  in  some  P7  as  defined 
in  3.3. 


Finally  let 


F7(x7)4F(x) 


Tj=r0l 


i(ty 


(3.4) 


where: 

(i)  F(x)  is  as  defined  in  equation  2.14  of  Chapter  1. 

(ii)  x  =  col(a.1,a2 . OLq.r) 

(iii)  x7  =  co^pu .02 . ^,.rTl.....rT  ) 

(iv)  71  denotes  the  i'th  element  of  the  fault  index  y 

(v)  r0i  is  the  nominal  value  of  the  i’th  parameter. 

If  a  system  fault  is  known  to  lie  in  P7  then  the  actual  values  of  the  faulty  parameters 
must  satisfy  F7(xy)  =  0.  Suppose  a  system  is  nf-fault  diagnosable,  then  there  will 
generally  be  one  fault  index,  y,  for  which  F7(x7)  =  0. 

4.  Limited  Fault  Algorithm 

This  section  has  the  following  objective: 

Given  a  circuit/system  which  is  nf-fault  diagnosable  and 
which  has  a  fault  characterized  by  a  parameter  vector, 
r  e  P7,  yef}( N,nf),  develop  a  solution  algorithm  which  will: 
i)  determine  which  circuit/system  parameters  are  faulty 
(find  7);  and  ii)  solve  for  the  values  of  the  faulty 
parameters,  i.e.  find  r  e  Py,  7en(N,nf). 

The  information  available  to  accomplish  this  objective  is  the  fact  that  the  parameter 

vector,  r,  and  the  associated  ambiguity  variables,  py,  i=  1,2 . q,  which  combine  to  form 

the  composite  unknown  vector,  x,  must  satisfy  F(x)  =  0  AND  the  constraint  that  reP 
where  P  =  o  Py.  The  latter  constraint  makes  it  impossible  to  utilize  the  Newton- 

7C  ll(S.Pf)  7 

Raphson  iteration  directly  To  illustrate  the  nature  of  the  problem  imposed  by  this 
constraint  consider  the  following  "brute  force"  approach  as  a.  solution  method: 
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Try  to  solve  for  Xy  in  F7(xr)=0  for  each  7efl(N,nf).  If  for  a 
particular  y,  a  solution  fails  to  exist,  y  denotes  an  incorrect 
fault  combination.  If  a  solution  does  exist  then  y  is  a 
potential  fault  combination.  If  only  one  y  yields  a  solution, 
the  fault  is  uniquely  identified. 

The  modified  Newton-Raphson  algorithm  developed  in  [1]  is  a  stable  method  for  testing 
each  F7.  The  number  of  such  equations  which  must  be  tested  is  CN;nf  where 

Cn^=  mi"(n-m)!  (4l) 

The  major  shortcoming  of  such  an  approach  is  that  as  N  gets  large  with  n,>2.  the 
number  of  possible  faults  becomes  extremely  large  resulting  in  excessively  long 
computation  time. 

The  excessive  number  of  fault  possibilities  for  large  systems  is  the  main  reason 
for  considering  the  alternative  approach  described  in  this  section.  For  clarity  we 
begin  with  a  brief  overview  of  the  idea.  Recall  that  the  available  information  for  a 
solution  algorithm  to  determine  the  fault  index,  y,  and  solve  for  the  parameter  vector, 
r,  is: 

F(x)  =  0  (4.2) 

and 


r  e  P  =  u  P_ 

TencNji,)  7 


(4.3) 


Suppose  we  view  equation  4.2  as  a  set  of  constraints  and  reformulate  the  restriction  in 
4.3  in  the  following  way: 


Define: 

<p(r.y)  &  Ir-r^  (4.4) 

where  (i)  reRN,  (ii)  r7eP7,  (iii)  7eCl(N,nf),  and  (iv)  |  •  ^  denotes  the  Euclidean  norm. 
Notice  that  since  equation  4.3  is  satisfied  if  we  require  that  p>(r,7)  be  a  minimum. 
In  other  words  <p,  the  objective  function  to  be  minimized,  is  the  distance  from  the 
solution,  r,  to  a  point  in  some  P7.  This  distance  is  zero  for  any  solution  that  can  be 
described  by  a  fault  index,  7en(N,nf).  Thus  the  problem  statement  becomes: 

minimize  <p{ r,y)  (4.5) 

r€R*  7en(N,nf) 

subject  to: 

F(x)  =  0  (4.6) 

To  develop  some  insight  into  the  solution  process  to  be  developed  consider  first  a 
solution  method  for  the  general  nonlinear  constrained  minimization  problem  given  by 
Luenberger[6j  and  known  as  the  "Gradient  Projection  Method".  This  method  proceeds 
in  the  following  manner: 
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(i)  Given  a  starting  point  use  some  solution  procedure  to  find  a 
"feasible  point"  that  is  a  point  which  satisfies  the  equality 
constraints. 

(ii)  Next  update  the  feasible  point  by  adding  to  it  a  component 
which  lies  in  the  space  tangent  (at  the  feasible  point)  to  the 
curve  defined  by  the  equality  constraints.  The  direction  of 
the  tangent  component  is  the  projection  of  the  gradient  of 
the  objective  function  into  the  tangent  space  (hence  the 
Gradient  Projection  Method).  The  updated  point  is  the 
value  along  the  tangent  projection  for  which  the  objective 
function  is  minimum, 

(iii)  Using  the  updated  point  as  the  starting  point  repeat  these 
two  steps  until  the  procedure  converges  to  a  solution  point. 

To  apply  this  type  of  procedure  to  the  solution  of  4.5  and  4.6,  modify  the  second  step 
of  the  Gradient  Projection  Method  in  the  following  manner: 

(i)  Let  ref>(N,nf).  For  each  -yeP  compute  the  tangent  space 

component  for  which  is  a  minimum. 

Note:  Although  this  requires  a  search  as  in  the  direct  approach 
the  time  required  can  be  kept  relatively  small  if  the 
number  of  elements  in  P  is  small. 

(ii)  To  choose  the  elements  of  T  consider  the  parameter  values 
at  the  feasible  point.  Any  parameters  which  are  close  to 
nominal  (say  within  10%)  consider  good.  This  will  eliminate 
some  fault  combinations  from  consideration  thus 
decreasing  the  number  of  elements  of  P  relative  to  the 
number  of  elements  of  fl(N,nf). 

(iii)  Another  criterion  which  will  limit  the  size  of  T  is  to  consider 
the  component  whose  parameter  value  at  the  feasible  point 
is  farthest  from  nominal  as  faulty. 

We  now  present  the  details  of  the  algorithm  to  solve  the  limited  fault  problem. 
This  algorithm  has  the  structure  shown  in  Figure  3.  The  first  step  of  the  algorithm  is 
to  choose  an  initial  guess.  The  best  available  information  about  the  solution  is  the 
nominal  values  for  the  circuit/system  under  test.  Therefore  let  x°  =  Xo- 

The  next  step  is  to  find  a  feasible  solution,  that  with  x0  as  an  initial  guess  find  a 
point,  xeRpq*N,  which  satisfies  F(x)  =  0.  Recall  that  the  modified  Newton-Raphson 
iteration  step  is[  1  ]: 

Jr(xk)dk  =  -F(xk)  (4.7) 

xkM  =  xk+Adk 

where  A  is  chosen  to  satisfy 

|F(xk+Adk)J|  <  |F(xk)|| 

JF(xk)  can  be  expected  to  have  a  nontrivial  null  space.  If  it  did  not  then  there  is 


Figure  3.  Basic  Algorithm  Structure 


sufficient  information  to  find  the  solution  without  the  limited  fault  assumption.  This 
means  that  4.7  does  not  have  a  unique  solution.  To  make  the  solution  to  4.7  unique 
requires  that  dk  be  the  "least  squares"  solution.  This  "least  squares"  solution  can 
theoretically  be  found  via  Moore-Penrose  pseudo-inverse[7]  but  in  practice  the 
solution  is  determined  using  software  such  as  1MSL  routine  LLBQF[8].  Implemented  in 
this  manner  the  modified  Newton-Raphson  iteration  provides  a  means  of  finding  a 
feasible  point. 

The  third  step  depicted  in  Figure  3  is  the  tangent  space  update  of  the  feasible 
point  to  minimize  the  objective  function.  To  assist  in  the  explanation  of  the 
implementation  of  this  step  define  the  following: 

Xfe=  the  feasible  solution  resulting  from  Step  2. 

rfe~:  the  parameter  part  of  xfe  (i.e.  the  last  N  entries  of  xIe). 

Vf(x;,.')4  the  matrix  whose  columns  span  Null[Jp(Xfe)]. 
vr(xre)^  the  matrix  consisting  of  the  last  N  rows  of  Vp(Xfe). 

Tp(xfe)=  the  tangent  space  of  the  surface  F(x)  =  0  at  the  point  Xfe. 
xt^  a  point  in  Tp(xfe). 
rt=  the  last  N  entries  of  xt. 

D7^  for  y  =  <it ,ia . inf>  this  is  the  identity  matrix  with  rows  i1(  i2 

etc.  deleted. 

Note:  D7  is  a  matrix  operator  for  the  parameter  vector  which  produces  a  vector 
consisting  of  the  entries  of  r  whose  indices  are  not  contained  in  y. 

The  objective  of  this  step  is  to  find  a  point  in  the  tangent  space  which  minimizes 
the  function  <p.  If  xt  is  a  point  in  Tp(xfe)  it  has  the  following  characterization: 

Xt  =  Xfe  +  VP(xfe)/S  (4.8) 

where  /9eRM  and  M=dim|Null[Jp(x,e)] j.  Thus  the  goal  for  this  step  becomes:  Find  0 
which  minimizes  <p.  Since  the  objective  function.  <p,  involves  the  parameter  alone  only 
the  last  N  equations  of  4.8  require  consideration.  Thus 

rt  =  rf,  +  Vr(xfe)/9  (4.9) 

characterizes  the  parameter  part  of  xt  used  in  the  objective  function.  Let  y  be  fixed. 
Note  that: 


min  <p(r,,y)  =  min  |r.-rJ» 

*t*TF(xf,rV  n  W**)1  7112 

=  min||D7(r,e+Vr(xfe)/3-r7)||2 
0eRM 


(4.10) 


But  since  D7r7  =  D7r0  this  become: 


min  C0(rtly)  =  minlA^-Bja 

xteTF(xf<1) 


(4.11) 


where 


Ay  =  D7Vr(>:f0)  and  B7  =  D7(r0-rfc) 


(4.12) 


For  a  given  and  y  the  matrix  Ay  and  the  vector  B7  are  known.  Thus  equation  4.1  ■  i:; 
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a  linear  least  squares  problem  readily  solved  via  techniques  discussed  previousIy[7,8]. 

Let  denote  the  vector  which  satisfies  4.11.  To  find  the  minimum  of  the 
objective  function  for  all  rt  and  y,  compute  for  each  y  and  select  the  one  which 
satisfies: 


min  lA^-Byb 

•yeOlN.iif) 


(4.13) 


As  discussed  earlier  performing  this  computation  for  every  possible  y  is  unnecessary. 
To  eliminate  some  fault  indices  from  consideration  construct  a  set  rcn(N,nf)  of 
candidate  faults  based  on  the  following  criterion: 

(i)  If  the  i’th  entry  of  rfe  is  sufficiently  close  to  the  nominal 
value  for  the  i'th  component  (e.g.  within  10%)  then 
eliminate  from  consideration  all  faults  containing 
component  i. 

(ii)  If  the  j’th  entry  of  rfe  has  the  greatest  deviation  from 
nominal  then  eliminate  from  consideration  all  faults  which 
do  not  contain  component  j. 

Instead  of  the  minimization  in  4.13  perform: 

mipjikyPy-Bjz  (4. 14) 

The  complete  algorithm  has  the  structure  shown  in  Figure  4. 


5.  limited  Fault  Algorithm  Examples 

The  purpose  of  this  section  is  to  present  two  example  problems  which  will 
illustrate  the  theory  and  limited  fault  algorithm  developed  in  this  chapter.  The  first 
example  is  based  on  the  12  parameter  circuit  of  Figure  5.  The  solution  algorithm  is 
employed  to  analyze  16  randomly  selected  faults  within  this  circuit.  The  second 
example  is  the  26  parameter  circuit  of  Figure  1.  Practical  implementation  of  the  fault 
diagnosis  procedure  for  this  circuit  requires  the  use  of  sparse  matrix  techniques. 
Although  the  development  of  sparse  matrix  algorithms  is  beyond  the  scope  of  this 
report  the  example  is  included  to  demonstrate  the  feasibility  of  the  algorithm  for 
larger  systems. 

Example  5.1:  Consider  the  circuit  of  Figure  5.  The  circuit  input  and  outputs  are: 

u,  =  V,  y,  =  ICl  ya  =  V0  (5.2) 

The  parameter  definitions,  nominal  values  and  component  transfer  functions  appear  in 
Table  13. 


Figure  4.  Detailed  Algorithm  Structure 
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Table  13. 

Component  information  for  Example  5. 1 


Parameter 

Nominal  Value 

Definition 

Zi(s,ri) 

<"i 

.1 

1  /c, 

ri/  s 

r2 

.5 

R* 

<*2 

<*3 

1. 

R„ 

r3 

<*4 

1. 

i/cM 

r4/  s 

r5 

.1 

l/c2 

r6/s 

re 

1. 

i/Rb 

r8 

<*? 

1. 

1/  Rg 

<"? 

<*a 

1. 

ma 

rBs 

<"g 

.1 

r8s 

<*io 

1. 

<*10 

<"11 

1. 

<"n 

<*12 

.5 

■BS 

r12 

The  nonzero  entries  of  the  sparse  set  of  connection  matrices,  Ltl,  Lig,  Lzi  and  Lzz 
appear  in  Table  14. 


Table  14. 

Nonzero  entries  of  the  connection  matrices  for  Example  5.1. 


Lit 

row.column 

value 

row.column 

value 

row.column 

value 

bn 

1.6 

1 

1,7 

1 

1.9 

1 

1,11 

1 

1,12 

1 

2,7 

1 

2,9 

1 

2,11 

1 

2,12 

1 

3,? 

1 

3.B 

-1 

3,9 

1 

3,10 

-1 

4,10 

1 

4,11 

1 

4,12 

1 

5,12 

1 

6,1 

-1 

7,1 

-1 

7,2 

-1 

7,3 

-1 

8,3 

1 

9,1 

-1 

9,2 

-1 

9,3 

-1 

10,3 

1 

11.1 

-1 

11,2 

-1 

1 1,4 

-1 

12,1 

-1 

12,2 

-1 

12,4 

-1 

12,5 

-1 

Lia 

6,1 

1 

7,1 

1 

9,1 

1 

11,1 

1 

12,1 

1 

Lai 

1.1 

-1 

1,2 

-1 

1,4 

-1 

1,5 

-1 

2,6 

1 

2,7 

1 

2,9 

1 

2,11 

1 

2,12 

1 

l'22 

1,1 

1 

The  nonzero  entries  of  I ,a,H  and  V,  computed  via  1MS1/B]  routines  LG1NF  and  LLSQF 
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respectively,  appear  in  Table  15. 


Table  15. 

Nonzero  entries  of  the  and  V  matrices  for  Example  5.1. 


matrix 

row,  column 

value 

row, column 

value 

l21r 

1.1 

-0.25 

2.1 

-0.25 

4,1 

-0.25 

5,1 

-0.25 

8,2 

0.2 

7,2 

0.2 

9,2 

0.2 

11,2 

0.2 

12,2 

0.2 

V 

1,1 

0.666025 

2,1 

-0.288675 

2,2 

-0.57735 

2,3 

-0.57735 

3,6 

1 

4,1 

-0.288675 

4,2 

0.768675 

4.3 

-0.211325 

5,1 

-0.288875 

5,2 

-0.21 1325 

5,3 

0.788875 

6.4 

0.861803 

8,5 

-0.138197 

8,6 

-0  138197 

6,7 

-0.138197 

7,4 

-0.138197 

7,5 

0.861803 

7.8 

-0.138197 

7,7 

-0.138197 

8,9 

1 

9,4 

-0.138197 

9,5 

-0.138197 

9,6 

-0.138197 

9.7 

0.861803 

10,10 

1 

11.4 

-0.138197 

11,5 

-0.138197 

11.8 

0.861803 

11,7 

-0.138197 

12,4 

-0.447214 

12,5 

-0.447214 

12,6 

-0.447214 

12,7 

-0.447214 

For  q=2,  Sj  =  jlO  and  sa  =  j.6  and  u(Sj)  -  u(sz)  =  1  we  compute  the  following 
nominal  otj: 


— .  157584e+0Q 

—  .353887e— 01 

— .33341 5e-t-00 

—  .130905e+00 

-,323898e+00 

—  .6274 19e-01 

.634977e  +  00 

-.  121971e— 01 

,30B373e-01 

+  j 

— .  105105e+00 

,803794e-01 

-,732278e— 01 

-,290531e+00 

.359565e+00 

.39B36Bc-01 

—  ,384998e— 01 

.384906e  +  00 

398368e+00 

39B368e-01 

-,3B4996e-01 

-,181031e+00  f-.167798e+00 

,399337e+00  -,395931e+00 

— ,247148e+00  -,?09381e-01 

.705033e+00  ,756794e-01 

,275746e+00  .  -.126582e+00 

^2  =  182022e+00  +  J  ,153923e+00 

— .  17B041e+00  -.116167e+00 

.213281e+00  -,420248e-01 

.252149e  —01  ,127968e+00 

.2132Ble+00]  -,420248e-01 


This  nominal  data  is  used  as  a  starting  point  for  the  solution  algorithm. 

The  next  step  in  the  example  is  to  simulate  several  3-faults  (3  parameters  differ 
from  nominal)  and  see  if  the  solution  algorithm  works.  Table  16  displays  a  set  of  16 
randomly  selected  fault  indices  and  corresponding  faulty  parameter  values.  Each  of 
these  faults  was  simulated  and  the  resulting  measurement  data  applied  to  the  solution 
algorithm.  The  final  column  in  Table  16  indicates  the  performance  of  the  algorithm  for 
each  fault.  In  11  of  the  16  cases  the  algorithm  correctly  identified  all  faulty 
parameters  and  determined  their  actual  values.  In  four  cases  two  of  the  three  faults 
were  correctly  identified  and  in  one  case  the  algorithm  selected  as  faulty  a  component 
which  was  actually  good. 

Table  16. 

Fault  list  and  algorithm  performance  summary  for  Sj=jl0.  and  S2=j.6. 


Actual  Fault  Index 

Faulty  Parameter  Values 

Algorithm  Fault  Index 

<6,7, 10> 

r6  =  .5 

r7=2. 

ra0=2. 

<6.7, 10> 

<2,5,9> 

r2  =  l. 

r5=.05 

r9=.2 

<2,9> 

<1,4,9> 

r,=  2 

r4=.5 

r9=.2 

<1,4,9> 

<2,5,6> 

r2=.3 

rs=.07 

re=1.4 

<2,6> 

<4,5,9> 

r4=2. 

rs=,05 

rB=.2 

<4,5,9> 

<2,6,1 1> 

r2=l. 

r8=.5 

rn  =  .5 

<2,6, 11  > 

<2,6, 12> 

r2  =  l. 

rs=,5 

r12=.25 

<2,6, 12> 

<7,8,9> 

r7  =  ,5 

rB=2. 

r0=.2 

<7,8,9> 

<2,4, 12> 

r2=.25 

r4=.5 

r12=.25 

<2,4, 12> 

<5.6, 8> 

r5=.0b 

re=.7 

re=1.6 

<6,8> 

<2,8, 10> 

r2=.25 

r0=2. 

rio=2. 

<2,8, 10> 

<3,4,7> 

r3=.5 

r4=2. 

r7=  2. 

<3,4,7> 

<3,6,9> 

r3  =  2. 

”1 

05 

II 

r9=  .2 

<3,6,9> 

<3,7, 12> 

r3  =  2. 

r7=2. 

ria  =  .25 

<3,7, 12> 

<2,4,5> 

r2  =  .3 

r4=1.5 

r5=.07 

<2,4> 

<5,7, 12> 

r5  =  .05 

r7=2. 

r  iz=  1- 

<7, 11,1 2> 

Although  the  circuit  is  theoretically  3-fault  diagnosable  the  results  in  Table  16 
indicate  that  there  are  some  practical  problems  associated  with  the  determination  of 
certain  faults.  For  any  given  fault  the  most  reliable  indication  of  the  theoretical 
capability  to  determine  the  fault  is  a  test  of  the  Jacobian  evaluated  at  the  fault  but 

fins  is  i : • . ; : ;  uciieol  since  il  is  impossible  In  anticipate  all  possible  faults,  fpstoad  the 
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diagnosability  test  is  based  on  of  the  rank  of  selected  columns  of  the  Jacobian,  Jp, 
evaluated  at  the  nominal  point.  Theoretically  the  results  at  the  nominal  point  hold  for 
almost  all  faults.  Unfortunately  it  is  quite  possible  for  a  matrix  which  is  theoretically 
full  rank  to  be  less  than  full  rank  for  a  solution  algorithm  due  to  the  finite  word  length 
of  the  computer. 

A  technique  for  circumventing  this  problem  is  to  perform  the  diagnosis  at  several 
sets  of  test  frequencies.  Only  those  faults  which  are  poorly  conditioned  at  all  test  I 

frequencies  used  would  not  be  detectable,  Suppose  that  the  present  example  is 
repeated  with  the  following  test  frequencies:  S!  =  j4.  and  s2  =  j.2.  The  resulting 
nominal  a*  are: 


0.1 


— .  152138e+00 
—  ,291921e+00 
-,309870e+00 
.639839c +  00 
.41 1483e-01 
,892647e— 01 
— .343296e+00 
,673079e-01 
,267029e+00 
.673079e— 01 


—  ,306290e— 01 

—  ,223071e+00 

—  ,464630e  —01 

—  ,987588e-02 
.  — ,825390e— 01 
J  ,328109e— 01 

,961794e— 01 

—  ,667572e— 01 
,269232e+00 

—  ,667572e— 01 


(5.5) 


,522551e— 01 

—  ,460309e+00 

,550190e+00 

•771991e— 01 

-.  196403e+00 

—  .B90725e— 01 

,577404e+00 

.307545e+00 

,304301e+00 

+  j 

,939940e— 01 

-.2461  62e+00 

—  ,347141e-01 

— ,898594e-01 

—  ,883199c— 01 

.201 163e+00 

.7BB47Be— 01 

-,157690c -01 

,402325e-01 

.201 183e+00 

,788478e— 01 

(5.6) 


The  results  of  the  diagnosis  of  the  same  16  faults  appears  in  Table  17. 
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Table  17. 

Fault  list  and  algorithm  performance  summary  for  st=j4.  and  sz=j.2. 


Actual  Fault  Index 

Faulty  Parameter  Values 

Algorithm  Fault  Index 

<6,7, 10> 

ra  =  .5 

rv=2. 

rio=2. 

<6,7, 10> 

<2,5,9> 

r2  =  l. 

rs=.05 

rB=  .2 

<2,5,9> 

:1,4,9> 

r,=  2 

r4=,5 

r9=.2 

<  1,4,9> 

<2.5, 8> 

rz  =  .3 

r5=.07 

r8=14 

<2,5,6> 

<4,5,9> 

r4=2. 

r5=.05 

r9=.2 

<4,5,9> 

<2,6, 11  > 

r2  =  l. 

re=5 

rn  =  .5 

<2,8,1 1> 

<2,8, 12> 

r2  =  l. 

re=.  5 

ri2=.25 

<2,6, 12> 

<7,8,9> 

r7  =  .5 

<r 

II 

lO 

r9=  2 

<5,7,9> 

<2,4, 12> 

r2=.25 

r*=  5 

r12=.25 

<2,4,12> 

<5,8,8> 

r0=.05 

re=.7 

r3= 1.6 

<5,6,B> 

<2,8, 10> 

r2=25 

ra=2. 

rio=2. 

<2,3, 10> 

<3,4,7> 

r3=5 

r4=2. 

r7=2. 

<3,4,7> 

<3,6,9> 

ra=2. 

r8=  5 

r9=  .2 

<3,6,9> 

<3,7, 12> 

r3=2. 

r7=2. 

r,2=.25 

<3,7, 12> 

<2,4,5> 

rz  =  .3 

r4=1.5 

r5=.07 

<2.4, 5> 

<5,7, 12> 

ra=.05 

r7  =  2. 

ri2=l 

<5,7, 11  > 

A  diagnosis  based  on  the  combination  of  the  results  of  Tables  15  and  16  identifies 
ALL  faulty  components  in  the  sixteen  randomly  selected  fault  combinations.  The  only 
diagnosis  errors  are  three  of  the  sixteen  cases  in  which  a  good  component  was 
identified  as  faulty 

Example  5.7:  Next  reconsider  the  circuit  of  Figure  1.  Due  to  the  large  size  of  this 
example  a  detailed  analysis  is  not  feasible  without  the  utilization  of  sparse  matrix 
techniques.  Although  the  adaptation  of  the  solution  program  to  sparse  matrices  is 
beyond  the  scope  of  the  current  research  a  limited  analysis  of  the  example  is  included 
to  demonstrate  the  feasibility  of  this  approach  for  large  circuits  and  to  illustrate  some 
of  the  properties  of  the  algorithm. 

The  choice  of  inputs  and  outputs  are  denoted  in  the  figure  as  Uj  and  yj 
respectively.  For  this  example  N  =  26,  n0=4  and  p=22.  Notice  that  five  of  the  seven 
accessible  nodes  are  utilized  to  provide  the  single  input  and  four  outputs  shown  in 
Figure  1.  Normally  it  would  be  desirable  to  utilize  all  accessible  points  to  acquire  the 
maximum  of  available  information  for  the  diagnosis  procedure  however  we  have 
elected  to  let  nf=5  and  wish  to  illustrate  the  case  where  n„<nf. 

The  nonzero  entries  of  the  sparse  set  of  connection  matrices,  L]2  and  Lgj,  appear 
in  Table  16  and  those  of  l,n  appear  in  Tabic  1  (Chapter  i).  All  entries  of  L22  are  zero. 
The  computation  of  L2,R  and  V  was  performed  via  the  IMSL[6]  routines  LGINF  and  LSVDF 
respectively.  These  too  are  spare  matrices  with  their  nonzero  entries  given  in  Table 

16. 


Table  IB. 

Nonzero  entries  of  matrices  for  Example  5.7. 


row, column  value  row,  column  value  I  row,  column  value 


2,9 

1 

3,19 

4,22 

1 

9,2 

1 

19,3 

26,3 

0.5 

3,3 

1 

4,4 

6,6 

1 

7,7 

10,9 

1 

11,10 

13,12 

1 

14,13 

16,15 

1 

17,16 

19,1 

1 

20,18 

23,20 

1 

24,21 

The  component  transfer  functions  are:  Zjfs.rj)  =  rjs  for  i=4, 5, 10, 11, 16, 17, 33,34  and 
Zi(s,rj)  =  n  for  the  remaining  i.  The  nominal  component  values  as  well  as  parameter 
values  for  three  5-faults  appear  in  Table  19.  Note  that  the  component  values  are 
scaled  to  improve  the  numerical  condition  of  the  problem.  The  impedance  scale  factor 
is  102  and  the  frequency  scale  factor  is  107. 


Table  19. 

Nominal  parameter  values  and  selected  5-faults. 


component 

units 

nominal 

— 

fault  1 

<4.13,18,20,22> 

fault  2 

<1,9,13,20,24> 

fault  3 

<3,8, 14,17, 25> 

i 

ohm 

12 

12 

16 

12 

2 

mho 

0.1 

0.1 

0.1 

0.1 

3 

ohm 

56,7 

56.7 

58.7 

35 

4 

farad 

50 

80 

50 

50 

5 

farad 

5 

5 

5 

5 

6 

mho 

10 

10 

10 

5 

7 

ohm 

30 

30 

30 

30 

8 

mho 

0.1 

0.1 

0.1 

0.1 

9 

ohm 

10 

10 

3 

10 

10 

farad 

50 

50 

50 

50 

11 

farad 

5 

5 

5 

5 

12 

mho 

10 

10 

10 

10 

13 

mho 

0.3 

0.2 

0.5 

0.3 

14 

ohm 

10 

10 

10 

15 

15 

ohm 

2 

2 

2 

2 

16 

farad 

50 

50 

50 

50 

17 

farad 

5 

5 

5 

8 

16 

mho 

10 

18 

10 

10 

19 

ohm 

10 

10 

10 

10 

20 

mho 

0.3 

0.5 

0.5 

0.3 

21 

ohm 

10 

10 

10 

10 

22 

ohm 

10 

8 

10 

10 

23 

farad 

50 

50 

50 

50 

24 

farad 

5 

5 

3 

5 

25 

mho 

10 

10 

10 

15 

26 

ohm 

0.7B 

0.78 

0.78 

0.78 

The  number  of  test  frequencies  is  q=2  and  the  inputs  and  test  frequencies  are 
u(sj)=u(j.2)=2.  u(Sj)=u(j.Ol)  =  .l  i 


The  nominal  values  for  the  Oj,  i=l,2  appear  in  Table  20. 


Table  20. 

Nominal  values  for  o^j  and  a? 


real  cart 


,658870c +00 
■956600e— 04 
.1609B8e+01 
.380b48e— 01 
.  182796e+01 
.956600e— 02 
,312299e+00 
.  103920e-03 
,571410e+00 
.  125561e+01 
.  103920e-01 
245900e+00 
.770900e— 02 
543437e+00 
.  l93858e+00 
.512714e— 01 
770900e-01 
.  1 98924 e +00 
.444493e— 02 
•254349e— 01 
.  173045e+00 
.444493e— 01 


,460866e+00 
.380548e— 03 
.272186e+01 
,956600e— 02 
,237069e+00 
,380546e-01 
.  12S581e+01 
,571410e-02 
.  103920e-01 
,312299e+00 
,571410e+00 
.  175069c +00 
,193858c -01 
,237413e+00 
,770900e-01 
,250221e+00 
.  193858e+00 
307628e+00 
.  25434 9e— 02 
,444493e— 01 
,351821e+00 
. 25434 9e -01 


real  Dart 


.  23607 2e +00 
.  155763e— 03 
. 91066 6e +00 
,329066c -03 
234965c -01 
.  155763e— 01 
,501926e+00 
,830079c -03 
,433474c -02 
■,530596e— 02 
,830079e— 01 
,434513e— 01 
.101418c -01 
.20755  le +00 
.134318c  -02 
,410939e-01 
.101418e+00 
,201038e-03 
.3571 17e— 02 
.  105925e— 03 
,452230e— 03 
.3571 17e— 01 


,357115e+00 
-,858132e— 04 
-,332734e+00 
,778815e— 03 
,270696e— 01 
— ,658132c-02 
.1061 19c+00 
,866948c -03 
•415040e— 02 
,250963c— 01 
.866948c -01 
,524036e— 01 
— ,268636e— 02 
—  .441227e— 01 
,507090e— 02 
,476014e— 01 
— .26B636e  —01 
,265215e— 01 
,21 1651e— 03 
.I78559e— 02 
,247148e— 01 
.21 1851e— 02 


The  results  of  algorithm  for  this  example  comprise  Table  21.  The  program 
correctly  identified  faults  1  and  2  and  four  of  five  faulty  components  in  fault  3.  The 
difficulty  in  identifying  parameter  14  is  another  example  of  a  poorly  conditioned 
problem  however  in  this  case  the  difficulty  is  not  frequency  dependent.  The  sensitivity 
of  the  four  outputs  to  changes  in  parameter  14  is  extremely  small  at  all  test 
frequencies.  Under  such  a  circumstance  it  is  reasonable  to  expect  difficulty  in 
determining  the  parameter  since  the  outputs  contain  little  information  about  it. 

Table  21. 

Results  for  example  5.7. 


fault  # 

actual  7 

1 

|  <4,13,16,20,22> 

solution  y 


<1,9,13,20,24>  <1,9, 13,20, 24> 


<3,6,14, 17,25> 


<3,6, 1 7,25> 


This  solution  is  a  good  example  of  the  reduction  in  the  number  of  fault 
combinations  which  results  from  the  determination  of  a  feasible  point.  For  this 
example  nf=5  and  N  =  26.  Consequently  there  are  C2os  or  65760  possible  fault 
combinations.  Thus  a  "brute  force"  approach  to  this  example  requires  the  solution  of 
657M0  different  sets  of  fault  diagnosis  equations.  One  the  other  hand  the  algorithm 
dr  vi  'i  p-  0  here  finds  a  finds  a  feasible  point  first  and  then  searches  for  the  fault 
among  the  parameters  winch  deviate  significantly  from  nominal  at  the  feasible  point. 
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For  example  in  the  solution  for  fault  1  there  were  10  parameters  which  deviated 
greater  than  ±10%  from  nominal.  If  the  remaining  components  are  considered  good 
and  the  parameter  with  the  greatest  deviation  from  nominal  is  considered  definitely 
bad  then  four  faulty  components  remain  to  be  found  from  nine  possibilities.  This 
means  that  the  algorithm  must  test  C94  or  126  different  faults.  This  is  significantly 
fewer  that  65780. 

These  examples  employed  a  FORTRAN  program  based  on  the  solution  algorithm 
described  in  the  previous  section  and  compiled  and  executed  on  a  VAX  11/780 
computer.  It  is  also  noteworthy  that  both  examples  identify  a  reasonable  number  of 
simultaneous  faults  while  achieving  the  goal,  n0<VT?;  set  in  [10]. 

6.  Conclusions 

Clearly  the  Tableau  Approach  is  readily  adaptable  to  the  assumption  of  a  limited 
number  of  simultaneous  faults.  The  major  problem  associated  with  this  assumption  is 
the  need  to  avoid  testing  the  enormous  number  of  fault  combinations  which  occurs  for 
large  systems  with  more  than  a  couple  of  faults  possible.  The  algorithm  developed 
here  avoids  this  problem  by  utilizing  the  information  available  at  the  surface 
described  by  the  equality  constraints.  The  most  significant  aspect  of  this  approach  is 
that  it  does  not  require  the  number  of  test  points  to  be  greater  than  the  number  of 
assumed  faults.  This  allows  a  reduction  in  test  point  requirements  over  other 
methods. 
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APPEND IX  A  Full  Oiasnosis  Pr o sram 


LATEST  VERSION - 6  DEC  32 - 

prosram  driver  <  i  nput.  output.  tape5=input.  tape6=out,Putl 


c  This  prosram  solves  the  fault  diasnosis  equations  under 

c  the  assumption  that  the  equations  are  quadratic  Sparse 

o  matrix  techniques  are  not  used.  Since  all  matrices  are 

c  s  par  s  e  •  a  s  par  s  e  ma  tr  i  x  i  m  p  1  em  e  n  ta  t  i  c  n  i  s  a  dv  i  sa  hie. 

o 

c  Development  stase  be  sun  21  April  1932. 

c  Finished  with  completion  of  Chap 5  on  .Julv  19.  1932. 

c  Modified  for  CAS  paper  October  S.  1982. 

c  Modified  to  solo e  for  normalised  parameters  December  o •  1‘ 

c  Modified  to  use  libs f  instead  of  llsqf  December  6 ■  1932 

i“: 

c  This  version  uses  the  Newton— Raphson  search  direction  but 

c  uses  the  fact  that  the  components  of  F(x)  are  quadratic  t 

c  i  n t er p o  1  a t e  MFC  x  >  I  I  **2  a  1  o n 3  the  s ear  c h  dire c t  i  o n.  S  i  n c e 

o  MF<  x )  M  ##2  i  s  f  o  ur  t  h  or  d er  a  f  i  v  e  p  o  i  n  t  i  n  t  er  p  o  1  a  t  i  o  n  i  s 


c 

ex  a  c  t. 

The  inter po 

lant  is  used  to  find 

t  h  e  p  o  i  n  t. 

a  1  o  n  j 

c 

o 

s  ear  c  hi 

d  i  r  e  c  t.  i  o  n  f 

or  wh  i  ch  MFC  x  )  1  1  **2 

i s  minimi 

red. 

o 

c 

c 

variabl 

e  list 

f: 

i  n  t  e  s  er  s 

i": 

i f: 

n 

number  of 

p a r am  e  t  er  s .  a  1  s  o  r  ow 

dime  n  s i  >  j  n 

o  f  Z 

c 

q 

n  urn  t<  er  o  t 

t  e  s  t  f  r  e  q  u  e  n  c  i  e  s 

f: 

p 

dime  n  s i o n 

o  f  t  he  am  hi  3  u  1 1 

i": 

C 

i  n 

d  i  m  e  n  s  i  o  n 

of  the  input  vector 

i.: 

ou 

dime  n  s i o n 

o  f  t  hi  e  o  u  t  F'  u  t  ,-..1  e  c  t  or 

i": 

S  p  OW 

p oilier  of  s 

in  Z  C  s  )  .  us  u-j  1  1  v  - 1  • 

0  ■  + 1 

i‘: 

o 

i  p 

inteser  w or k  vector  used  bv  ims 

1  r  o  u  1 1  n  e 

1 1  s  q  f 

i*: 

kb 

u  s  e  d  b  .•  i  m 

s 1  r o  u  t i n  e  11 s q f 

r  ea  1 
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ccm  connection  matrix 


com  connection  matrix 


ccm  connection  matrix 


ccm  connection  matrix 


121r 


r i 3 h t  inverse  of  121 


matrix  of  basis  vectors  for  the  null  space  of  121 
nominal  parameter  values 


r  a  c  t 


actual  parameter  values 


rkth 


parameter  values  during  the  iteration  process 


rkthnm  rkth/rnom  (normalised  unknown) 


deviation  from  nominal  for  actual 


real  matrix  equivalent  to  if 


equivalent  to  f  (see  f) 


equivalent  to  d  (see  d) 


equivalent  to  x  (see  x) 


used  by  imsl  routine  llsqf 


n  or m s  q  function  nam  e 


lambda  fun c t i o n  name 


array  of  function  va lues  along  the  search  d i r e c t i o n 


timel 


beginning  c  p  u  time 


t  i  m  e2 


ending  c p  u  time 


c  om  p  1  ex 


array  of  test  frequencies  (equiv  to  ss) 
array  of  input  vectors  ( equiv  to  uu) 
right  hand  vector  of  algor  (equiv  to  ff) 
search  direction  (equiv.  to  dd) 
point  along  search  direction  (equiv  to  xx) 


temporary  storage 


work  storage 
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c  Yinnom  array  of  nominal  test  outputs  (equiv  to  yymnom) 

c 

c  alpha  array  of  nominal  ambisuity  vectors  (equiv.  to  aalpha) 

c 

c  ymact  array  of  actual  test  outputs  < equiv.  to  yymact) 

c 

c 

o  Since  the  free-format  read(5,  *)  won't  read  complex  numbers 

c  the  variables  s,  u»  ymnom,  al  pha.  ymact  are  made  equivalent  to 

c  real  arrays  of  double  si2e.  The  reals  are  read  but  the 

c  complex  are  used, 

c 
c 
c 

logical  convrs 

inteser  n,  q,  p,  in<  ou>  2 pom.' (30) 

inteser  Ida,  ld.if,  ldlll,  1  daO,  ldbO,  ldal,  ldv,  ip(200)>  kb 
real  111(30, 30), 112(30, 5), 121(10, 30), 122( 10, 5). 121r(30. 10) 
real  v<30, 25), rnom(30) , ract(30) , a(250, 250),  ff(250),  dd(250) 
real  ss< 10) ,  uu( 10, 5) , vvmnom(20, 5) >  aal Pha( 50,  5) . yymact < 20,  5) 
real  normsq,  lambda,  tol,xx(250),  3(5),  dev(30) 
real  timel, time2, r kth< 30) , rkthnm( 30) ,  c(4) 
complex  s( 5) , u( 5. 5) , ymnom( 10, 5) , al Pha ( 25, 5) , ymact ( 10. 5) 
complex  ..if  ( 150,  250),  d(125),  f(125),  b0<30,  5) ,  a0(  30,  5) 
compl ex  tmp(200) , aa(30,  5) ,  x( 125) 
equivalence  (  ss,  s) ,  (  uu,  u) ,  (aal  pha,  al  pha) 
equivalence  ( y yma c t , Yma c t ) •  (yymnom, vmnom) 
equivalence  (  f  f ,  f ) .  <  dd,  d) ,  ( xx,  x) 
c 
c 

call  second( timel ) 

c 

c  set  the  row  dimensions 

c 

1 da=250 
1 dJ  f=150 
1  dll  1=30 
1 dv=30 
Ida 1=25 
1 da0=30 
1 db0=30 

o 

c  set  precision  of  solution  i  e  .01  =  1 V. 

prec  =  .  001 

c 

c  I nput  secti  ■  r 

c  read  in  from  standard  input 

c 

read(5,  *)  n,  q,  p.  in,  ou 
read<5,  *)  (rn om(i),i  =  l,n) 
r  ea  d ( 5 • * )  <  ra  c t ( i ) ,  i  =  1 ,  n ) 
do  55  i = 1 , n 

55  r ead ( 5,  * )  <  1 1 1  (  i ,  j  ) ,  j  =  1 ,  n ) 

do  60  j  =  1 ,  i n 

r  ea d ( 5,  * )  <  11 2  <  i ,  i ) ,  i  =  1 .  n ) 
do  65  i=l,  ou 


60 
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65 

70 

75 


80 


85 

90 

95 

100 


98 


o 

e 

c 

101 

102 

103 

107 

no 


111 

112 

113 


r  ead<  5.  *)  (  121  (  i,  J  >  >  i=l,  n) 
do  70  J=1 . i n 

read  (5.  *) < 122< i>  • ) . i  =  1 » ou) 
do  75  ,i  =  l»  ou 

read<5.  *)ei21r<i.  .i ) .  i  =  1 .  n) 
do  80  j *  1 .  p 

read<5.  *)  <m<  i,  .<) .  i  =  l.  n) 
read<5.  *)  (  2 p om.i (  i  )  .  i  =  l*  n) 
r ea d ( 5.  * )  <  s  s  e  i ) .  i  =  1 .  2*  q ) 
do  85  i=l . i n 

readeS.  *)  <  uue2*i-l.  i )  .  uu(2*i»  .• ) .  -i  =  l.  q) 
do  90  i  =  l*  ou 

read<5# * ) < YYmnom< 2*i-l .  ,• ) »  YYmnome2*i.  i ).  ,i  =  l»  q) 
do  95  i  =  l .  p 

read< 5.  * ) (aal phae 2*i-l .  J  ) .  aal  pha< 2*i »  J  ) .  i  =  l .  q) 
do  100  i=l .  ou 

read< 5.  *) ( YYmacte 2*i-l .  ,i ) ,  YYmacte  2*i<  i )  >  ,i  =  l.  q) 

nra  p  h=0 
i  tmax=200 
do  98  i=l»  p 
r eade 5-  *.  end=101 ) 
r ead< 5.  *.  end=101 )  nraph.  itmax 


print  a  program  heading 
write  <6.  102) 

f orma t <  //5x .  38 hFa u  1 1  Diagnosis  Pro gram-T a b  1  ea u  Method*/. 
lOx.  28h<  Interpolate  to  find  lambda)) 
wr  i  te<  6.  103) 

f  orma t  < / 1 Ox. 18 hN  om i na 1  Par am  e  t  er  s ) 
wri  te<6.  107)  <  "r<  ".  i,  y  >  =  '.  rnom(  i  ).  i  =  l.  n) 
formate 4<3x.  a2.  i2.  a2.  elO  4)  ) 
wr i t  e ( 6 .  110) 

formate /18x. 16hTest  Frequencies) 

do  111  i  =  1 .  q 

wri te<  6. 1 12)  i,  se i ) 

wr  i  t  e  ( 6 .  1 1 3 )  <  u  u  ( ,i ,  i  ) .  i  =  1 ,  2*  i  n ) 

c  o  n  t  i  n  u  e 

formate  12x.  2hse,  i2.  2h)=.  ell  4.  3h  +,«.  ell  4) 
formate  12x.  2hu=.  10 e  f4  1,  lx.  f4  1 . 5x)  ) 


c  Compute  the  Meet  or  bO  =  1 2 1  r  £vma  c  t- 1 22*  u> 

o 

do  130  i  =  1  >  q 
do  1 25  j  —  1 .  o u 
tm p  <  i )  =  cm  p  1  x  e  0  •  0.  ) 
do  120  k=l. in 

120  tmpe  ,i  )=tmp(  j  )  +  122e  j.  k)*u<  k.  i ) 

1 25  tm p  e  i )  =  Yma  c  t  <  .i ,  i )  -  tm p  e  j ) 

do  130  j  =  1  •  n 
bO  e  i ,  i  )  =  cm  p  1  x  e  0  .0  ) 
do  130  k=l.  ou 

130  bOe  i.  i )  =  bO  (  ,i.  i  )  +  121re  ,i .  k)*tmp<  k) 
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c  Compute  the  vector  aO 

c 

do  204  i  =  l  > e 

do  204  j  =  l,  n 

aO(  ,i.  i  )  =cmplx<0.  .0  ) 

do  201  k=l, n 

201  aO<  j,  i )=aO<  i.  i  )  +  l  1 1  <  j,  k)*bO<  k,  i  ) 

do  204  k=l » in 

204  aO  <  j ,  i  )  =aO  <  ..i ,  i  )  + 1 1 2  <  j ,  k )  *  u  (  k,  i  ) 

c 

c 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

do  280  i  =  l,  n 
r  k  t  h  nm  <  i  )  =  1 . 

280  dev< i )=<ract< i )-rnom< i ) )*100.  /rnom( i ) 

i ter=0 
xlam=-l. 

290  i ter-i ter+1 

if  (iter  at.  itmax)  so  to  900 

c 

d  o  295  i  =  1  >  n 

295  rkth< i ) =r kthnm< i ) *rnom< i ) 

c  eve  1  ua  t  e  f 

do  300  i=l > e 

call  e  ua  d  f  <  f  <  i  *  n-  n+ 1 ) .  r  k.  t  h .  a  1  p  ha  <  1 »  i ) ,  s  (  i  ) ,  zpow,  111.  1  d  1 1 1 , 

S<  v,  ldv,  aO(l>  i).  bO(l>  i),  n,  p,  tmp) 

300  continue 

i  f  (  normss<  f  f .  2*n#s ) .  le.  1.  e-1 0)  so  to  900 

o 

c 

c  stop  because  the  algorithm  appears  stuck  in  a  relative  min. 

if(xlam.  ee.  0  )  wr  i  te<  6-  305) 

305  f ormat< lOx, 37hcausht  in  a  point  of  relative  minimum) 

if(xlam.  es.  0.  )  stop 

c 

cal  1  jacob(  .i  f,  1  d  i  f,  rkth.  al  pha,  1  dal ,  s,  n>  p>  e,  spowj  111.  1 dl 1 1 , 

&  v.  1  dv,  aO.  IdaO.  bO.  ldbO.  tmp.  aa) 

c 

call  cxtorl(a>  Ida.  .if.  ldjf.  n.  p.  e) 
c 

cpar  since  Jacob  computes  C  partial  f  /  partial  rkth  1 

c  the  last  n  columns  of  the  ma  tr i x  a  must  be  a d i  u s t e d  to  s e t 

c  C  Partial  f  f  partial  rkthnm  3 

do  315  i  =  1 . 2*n* s 
do  315  i  =  l.  n 

315  a  (  i ,  2*  p*  s+  J  )  =a  <  i  .  2*  p*  s+  ,i )  *r  n om  <  i ) 

cpar 

c  see  imsldoc  lls^f  for  meaning  of  tol 

tol=0 

c  value  for  kb  indicates  a  is  assumed  to  have  independent,  cols. 

kb=2*p*e+n 
kb=0 

c  use  llsef  to  find  the  search  direction 

c  llsef  is  from  the  imsl  library 

cxxx  call  llsef(a.  Ida,  2*n*e,  2*p*g+n,  f,  tol.  kb,  d.  Imp,  ip.  ier) 
cddd 


nf=l 
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cddd 


c 


o 


C 


319 

321 

c 

C 


329 

331 

c 

339 


c 


400 


410 


420 


456 


it 


i  nd*l 
c  <  1 )  =0. 
c<  2)*0. 

C<3)  =  1. 

call  1 1  bsf  < a.  1  da.  '2*n*s,  2*p*s+n.  f >  Ida.  nf .  i  nd.  c.  d.  Ida.  ip.  ■* f  •  i  er ) 
NOTE:  the  r  part  of  d  is  normalized 

option  to  skip  the  lambda  determination  and  set  it  to  -1 

if  (iter.  le.  nraph)  so  to  400 

do  319  J  =  1 .  s 
do  319  i*l. p 

x< ( J-l )*P+i )=alPha( i.  j)-0.  0*d< < J-1 >*p+i ) 
do  321  i  =  1 »  n 

xx ( 2*p*s+i )  =  <  r kthnm( i ) -0.  0*dd<  2*p*  e+i ) ) *r nom( i ) 
the  —0.  0  is  the  lower  end  of  the  interpolation  data 
compute  function  If  1**2  at  lambda=0  -  25  5  -  75  -1 

do  339  ilam=l>5 
d  o  333  i  =  1 .  s 

call  suadf< f ( i*n-n+l ) .  xx < 2* p* 1 ),  x< ( i-1 )*p+1 ).  s(  i ).  zpow. 

111.  1  dl  1 1 ,  u.  1  du.  aO(  1 .  i  ) ,  bO<  1 ,  i  ) .  n,  p.  tmp ) 

c  o  n  t  i  n  u  e 

s( i lam) =normss< f f , 2*n*s) 
do  329  J  =  l,  s 
do  329  i  =  l . p 

x ( ( j-1 )*p+i )=x< < j-1 )*p+i )-  25* d( ( i-1 )*p+i ) 
d  o  33 1  i  =  1 .  n 

xx <  2* p* s+ i ) =xx( 2*p*s+i )-.  25*dd< 2*p*s+i )*rnom< i ) 
the  -.25  is  the  increment  between  the  interpolation  data 
coriti  nue 
xlam=lambda< s) 

update  the  alpha's  and  r 


continue 
do  410  j  =  l.  i 
do  410  i = 1 .  p 

a  1  p  ha  (  i .  j )  =a  1  p  ha  <  i .  j )  +x  1  am*  d  (  (  —  1 )  *  p+  i  ) 

Conors*  true, 
do  420  i=l.  n 

i  f (abs(xlam*dd(2*p*s+i  ) )  st.  prec)  convrs*.  false. 
rkthnm(  i  )=rkthnm(  i  )  +xlam*dd<  2*p*s+i ) 
if(convrs)  so  to  900 

write  the  parameters  as  a  prosress  report 
write <6. *>  '  rank  of  a  :  int(c(4) ) 

i f  < i t  er .  s  t  nr a  p  h ) then 

wr i t  e  <  6 . 455 )  iter. x 1 am 


el  se 

wr i t e  <  6. 456 )  iter. x lam 

endi  f 

f  or ma t ( /5x . 26  h  n  or ma 1 i z  e  d  r  at  it  era t i o n  .  i 3 • 4 x . 
Sh<  lambda*.  el2.  5.  1  h )  ) 

f  or  ma  t  < /5x , 26  h n  orma 1 i z  e  d  r  at  iteration  .  i 3. 4x • 
8h<  lambda*.  el2  5.  Sh:  forced)  ) 
wr i t  e  <  6. 457 )  ( r k  t hnm  < i ) •  i  =  1 •  n ) 
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457  format(5x,  4el5.  5) 

so  to  290 
o 

c  write  out  final  values 

900  continue 

call  second<  time'2 ) 

wri te< 6, 923)  iter-1,  time2-timel 

923  format(/5x,  i5,  2x,  16h  iterations  and  »  el2.  6,  14h  sec.  required) 

wri te<6,  925) 

925  f ormat< /5x,  36hParameter  Values  at  Solution  Point  > 

Sc  19h%  dev.  from  nominal,  17h  V.  error  in  sol.  ) 

do  950  i=l,  n 

927  formats  lOx,  2hr<,  i2,  4h)  =  ,  ell.  4,  18x,  f  10.  2,  5x,  f  10  2) 

950  write<6,927)  i, rkthnm< i )*rnom( i ) . dev< i ) , 

S<  100.  *<  ract(  i  ) -r kthnm<  i  )  *r nom<  i  ) ) /ra c t <  i  ) 

stop 
end 


c 

c 

o 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  euadf<  f,  r,  alph,  s,  zpow,  111,  ldlll,  v,  1  dv,  aO,  bO,  n,  p.  tmp) 

This  routine  computes  the  nonlinear  function  which  corresponds 
to  the  fault  diagnosis  equation  at  one  test  frequency. 

argument  list: 


f 

output: 

f  =  Z  <  s  , r)CL  Valph  +aO< s  )  3  -  Valph  -bO< s  ) 
i  1 1  i  i  i  i 

r 

i  nput: 

par  am  e  t  er  v  e  c  t  cr¬ 

al  ph 

i  nput: 

am  b  i  3 u i tv  v e  c  t  or 

s 

i  nput: 

test  frequency 

2POW 

input: 

exponent  of  s  in  Z 

111 

i  nput: 

com  connection  matrix 

ldlll 

i  nput: 

row  dimension  of  111  in  ca 1 1 i n s  pr o sr am 

V 

input: 

null  space  basis  of  121 

1  dv 

i  nput: 

r ow  dimension  of  v  in  calling  pr o sr am 

aO 

i  nput: 

particular  solution  for  a 

bO 

input: 

particular  solution  for  b 

n 

p 

tmp 


input:  number  of  parameters  St  row  dimension  of  Z 

input:  number  of  columns  in  u 

input:  temporary  storage  <at  least  2#n) 


i  nteser  zpow<  1 ) .  n,  p,  1 dl 1 1 »  1  dM.  i ,  J 

real  r < 1 ) .  1 1 1 < 1 dl 1 1 . 1 ) , v< 1 dM. 1 ) 

compl  ex  al  Ph(  1 ) .  s,  f < 1 ) , aO( 1 ) .  b0< 1 ) .  tmp<  1 ) 


do  20  i  =  l.  n 
tmp<  i )  =0. 
do  10  j  =  l»  p 

tmp < i >  = tmp < i ) +m< i . J ) *al ph <  J ) 
f  < i )=-tmP< i >-b0< i ) 
do  40  i  =  l.  n 
tmp(  i  +  n)=0. 
do  30  ,i  =  l<  n 

tmp( i  +  n)  =  tmp< i  +  n)  +  lll ( i,  j )*tmp<  j ) 
tmp( i+n)=tmp< i+n)+aO< i) 
f ( i )=f ( i )+r< i )*s**2pow( i )*tmp< i+n) 
r  e  t  ur  n 
end 


subr  outi  ne  iacob<  •  f,  1  dJ  f »  r .  alpha.  1  dal .  s.  n.  p.  zpow.  111.  1  dl  1 1 . 

v.  ldM.  aO.  IdaO.  bO.  ldbO.  tmp.  aa> 


This  routine  computes  the  fault  diagnosis  Jacobian  to  be 
used  in  the  Newton-Raphson  or  other  iteration  schemes. 


ar  a  urn  e  n  t  list 

jf  output:  the  Jacobian  matrix 

ldJf  input:  the  row  dimension  of  Jf  in  the  calling  program 


r 

alpha 
1  dal 


i n p u  t :  pa r am  e  t  e r  v e  c t  or 

input:  ma tr i x  of  am b i 3 u i tv  oect or s 

i  n p u t  r ow  dimension  of  alpha  i n  t he  calling  pr o gram 


input:  array  of  test  f r eguencies 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C! 

0 

c 

c 

c 

c 

c 

c 

c 

c 

c 

0 

c 

c 

c 

c 

c 

c 

c 


c 

C 

c 

c 

c 


301 


304 

c 

c 

C 

c 


n 

P 

q 

2POW 
111 
1  d  1 1 1 

M 

1  dv 

aO 

IdaO 

bO 

ldbO 

Imp 

aa 


input:  number  of  parameters  &  row  dimension  of  Z 

input:  number  of  columns  in  v 

input:  number  of  test  frequencies 

input:  exponent  of  s  in  Z 

input:  ccm  connection  matrix 

input:  row  dimension  of  111  in  call  ins  prosram 

input:  null  space  basis  of  121 

input:  row  dimension  of  m  in  call  ins  prosram 

input:  array  of  Particular  solutions  for  a 

input:  row  dimension  of  aO  in  the  call  ins  prosram 

input:  array  of  Particular  solutions  for  b 

input:  row  dimension  of  bO  in  the  call  ins  prosram 

input:  temporary  st erase  (at  least  2*n) 

input:  work  array  (at  least  n*q) 


Note:  For  arrays  al  pha(  i »  i ) ,  a0(  i ,  ,i ) ,  b0<  i ,  ..i )  &  aa(i,-.i>  the 
subscript  j  means  the  data  of  this  column  corresponds  to 
test  frequency  s(i). 

inteser  2Pow(ldaO)»  n>  p,  q,  ld,if#  Idal,  ldlll,  1  du,  IdaO,  ldbO 
i  nteser  i»  ■>,  k,  J  .j 

real  r(ldaO),  111(1  dill,  1  dl  1 1 ) .  v(  1  dv, ldal) 

complex  ,i  f  <  1  d  J  f ,  1  d-.i  f ) ,  al  pha(  1  dal ,  q) .  s<  q) ,  a0(  1  daO,  q>  ■  b0<  1  dbO.  qi 
complex  tmp(  1  d-i f )  •  aa(  n.  q) ,  S2 


Compute  the  vector  aa-Ll l*V*al Pha+a0 

do  304  i*l ,  q 
do  301  ■>- 1,  n 
tmp(  i  >  =  cmplx(0.  ,0  ) 
do  301  k«l. p 

tmp(  i  )  =  tmp(  <  )+w(  f,  k)*al  pha(  k,  i  ) 
do  304  j  =  l  ,  n 
aa<  j,  i  )  =a0(  ■>>  i  ) 
do  304  k=l,  n 

aa(  i,  i  >=*aa<  i,  i )  +1 1 1  (  J,  k) *tmp(  k) 


Compute  the  Jacobian  JF(x> 


do  601  i*l ,  n*q 


do  601  J  =  l>  q*p+n 
J  f  (  i>  J  )  =  cmplx(0.  ,  0.  ) 
do  610  k=l»  q 
do  605  i = 1 > n 

Compute  dfi/dr  *  si(alphai) 

S2=s<  k)**2Pow< i ) 

Jf(n*<k-l)+i>  q* p+ i ) =aa < i ,  k )  *  s  2 
do  605  J  =  l,  p 
d  o  605  j  J = 1  >  n 

Compute  2 ( si . r ) *L1 1*V-V 

..if  <  n*<  k-1 ) +i »  p*(  k-1 )+..»)  = 

J  f  (  n*(  k-1  )  +  i»  p*<k-l)+J)+lll<i>  J  J  )*v<  J  J,  J ) 
do  610  i  =  l,  n 
do  610  j  =  l»  p 
s  2=  s  (  k )  **  2  p  ow  (  i  ) 

..if  (  n*<  k-1  )  +  i,  p*(  k-1 )  +  J  )=  ..if  (  n*(  k-1  )  +  i»  p*(  k-1  )  +  .i  )*r<  i  )*s2-v(  i.  J ) 


return 

end 


subroutine  cxtorl(a>  ldai  if.  ldif.  n>  p.  q) 

This  routine  converts  the  complex  Jacobian  into  an 
equivalent  real  matrix.  Note:  The  conversion  is  somewhat 
of  a  hybrid  since  the  original  unknown  vector  is  part 
complex  (the  alpha's)  and  par t  real  (the  r y s ) 


argument  list 


a 

0  u  t  put: 

the 

1  da 

i  nput: 

the 

Jf 

input: 

the 

Id. if 

i  nput: 

the 

n 

i  nput: 

the 

p 

i  nput: 

the 

q 


real  equivalent  of  Jf 

row  dimension  of  a  in  the  call  ins  prosram 
complex  matrix  to  be  converted 
row  dimension  of  Jf  in  the  call  ins  prosram 
dimension  of  r*  the  real  sub— vector- 
dimension  of  each  complex  sub— vector 


input:  the  number  of  complex  sub-vectors 
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integer  1  da*  ldJf,  n»  p.  a,  i,  i 
real  a. (  1  da*  1 ) 
compl  ex  J  f  <  1  dJ  f,  1 ) 
c 
c 

do  820  i  =  l . n*g 
do  810  j=l,  p*e 

a<2*i-l,  2*  J-l )=real <  jf  <  i,  J ) ) 
a  <  2*  i ,  2*  j  )  =r  ea  1  (  j  f  (  i ,  j  )  ) 
a  <  2* i - 1 ,  2* i ) =-a i ma  a  <  j  f  < i »  J ) ) 

810  a< 2*i ,  2* j-1 ) =aimag< j f < i * J ) ) 
do  820  ,i  =  l»  n 

a  ( 2*  i  - 1  ,  2*  p*  ■*+  j )  =r  ea  1  <  j  f  <  i  >  p*  g+  j ) ) 
820  a<2*i>  2*p*e+-j  )=aimas<  -i  f  <  i,  p*q+..i )  ) 

c 

r  e  t  ur  n 
end 


c 

c 

c 

real  function  normsi(  f.  n) 

c 

c 

inteaer  n,  i 
real  f  <  n ) 

o 

c 

n  or  m  s  e=0. 
do  10  i  =  l,  n 

10  normse»normsa+f < i )*f < i ) 
ret  ur  n 
end 

c 

c 


c 

o 


c 


c 


c 

c 

c 

c 

c 

o 


real  function  lambda(s) 

This  program  uses  5  points  along  a  function  g(  )  to  find  the 
point,  lambda,  at  which  the  function  is  a  minimum  The  ordinal, 
points  are  assumed  to  be  0  0.  25,  5,  -  75  and  -1 

If  S<x)=al*x**4+a2*x**3+a3*x**2+a4*x+a5  then 
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10  0  0  0  II  1*11  lj(-0.0)l 

13.  90625 e-3  -1.  5625e-2  .  0625  -  25  II  Ia2l  |3<-  25)  I 

I  .0625  -.125  .25  -.5  11x1*31  =  !  s(-0.  5)  I 

1.31640625  -.421875  .5625  -  75  11  1*41  l9<-75)l 


-1  II  1*51 


I  3<  — 1.  0)  I 


variable  list 


r  eal 


vand  *  *  =  s 


input:  contains  s<-0.  0) . s<-l.  0) 


va  n  d 


d*  t* : 


contains  the  Vandermonde  matrix 


vand2  work:  scratch  copy  of  vand  to  compute  a 

note:  It  would  be  more  efficient  to  use  the  inverse  of 
vand  here. 

a  work:  coefficients  of  the  interpolant 

becomes  the  coefficients  of  1st  derivative 

complex 


work: 


the  zeroes  of  the  1st  derivative 


subroutines 


leetlf  linear  equation  solver  from  imsl 


zpolr 


polynomial  root  finder- 


in  t  e  s  er  i ,  j 

real  vand (5*  5) , vand2<  5, 5) ,  a < 5) ,  3(5),  wk< 10) 
complex  2(3) 

data  van  d/0.  ,  3.  90625  e-3, .  0625,  .  31640625,  1  .  0  »  -1.  5625e-2,  -  125. 

-.  421375,  -1.  ,0  ,  0625,  25,  5625,  1.  ,  0.  •  -.  25,  -.  5,  -  75,  -1.  . 

1.  ,  1.  ,  1.  ,  1.  ,  1.  / 

set  fmin  to  the  value  of  the  function  at  the  current  point  so  that 

onlv  points  representing  a  decrease  from  this  value  can  be  selected 

fmi n=s< 1 ) 

do  80  i  =  l,  5 

do  80  j=l , 5 

vand2(  i,  ..i  )=vand<  i,  j  ) 

call  leatlf  ( vand2,  1,5,5,  s,  0,  wk.  ier) 


print  the  coefficients 
wri te<6,  100) (5-i,  s<  i ) ,  i  =  l>  5) 
format<//5<  lha,  i  1 .  1  h=,  elO  4,  2x)) 
ir oot=0 

find  coef.  of  1st  derivative 

do  10  i=l, 5 

a( i )=f 1 oat<5-i )*s< i) 

find  2*ros  of  1st  derivative 

call  zpolr < a.  3,  2,  i er ) 


do  50  i  =  l,  3 

i  f  <  cabs<  2<  i  )  ) .  at.  1.  el5)  so  to  30 
i  f  <  abs<  aimas<  2<  i  )  ) ) .  at.  1.  e-5)  30  to  30 

this  prevents  *  search  in  the  opposite  direction 
i  f  (real  <  2<  i  ) )  at.  0  )  so  to  30 

f lam=real <  2< i ) ) 
fma  9=0. 
do  25  k=l>  5 

fma  9=  fma  9+  9  <  k )  *  f  1  am**  <  5-  k  > 

if(fmin.  1  e.  fmas)  so  to  30 

iroot=i 

fmi n=fmas 

continue 

continue 

if(iroot.  ne.  0)  so  to  90 
i  f  <  9<  1 )  It.  0  )  then 

if  the  interpolant  is  upside  down  trv  lambda=l 
lambda=-l. 
ret  ur  n 

el  se 

lambda=0. 
wri te<6. 99) 
endi  f 

f ormat< 5x. 20hERR0R:  No  Minimum  !!) 
lambda=real <  2< iroot) ) 
if  <abs(  lambda),  at.  15.  )  then 
lambda=-l. 

write<&»  *)  'error:  lambda  too  large-' 

endi  f 
ret  ur  n 
end 


subroutine  second* time) 
real  time 
time=cputim< ) 
return 


APPENDIX  B:  Limited  Feu It  Diagnosis  Program 


EXPERIMENTAL  VERSION  - 11  DEC  82 

pr  o  am  n  f  ma  i  n 


This  program  solves  the  fault  diagnosis  equations  under 
the  assumption  that  the  equations  are  quadratic  and  that 
the  number  of  faults  is  limited. 

Development  stage  begun  December  6,  1982. 

Based  of  the  fault  diagnosis  program  driver,  f. 

logical  Conors 

integer  n.  q,  p.  in.  ou.  2 p ow ( 30 ) 

integer  Ida.  ld-jf,  ldlll.  1  daO.  ldbO.  ldal,  ldv.  ip<200) 

real  111(30,30).  112(30,5),  121(10,30),  122(10,5),  121r(30,  10) 

real  v(30, 25) . rnom<30) ,  ract<30) • a<250, 250) ,  ff(250),  dd(250) 

real  ss( 10) ,  uu( 10,  5) , YYmnom(20» 5) , aal pha(50. 5) . rvmact(20, 5) 

real  normss,  lambda,  xx(250),  g(5), dev(30). beta(30) 

real  timel, time2, rkth( 30) , rkthnm(30) , c( 4) . anul 1 ( 250, 360) 

complex.  s(  5) .  u<  5,  5) .  vmnom(  10.  5)  •  al  pha(  25,  5 ) .  vmact(  10,  5) 

complex  if < 150, 300), d(125),  f(125), b0(30. 5) , a0(30. 5) 

c  om  p 1  ex  tm  p ( 200 ) . aa ( 30 , 5 ) . x ( 125) 

equivalence  <  ss,  s) ,  (  uu»  u) ,  ( aal  pha,  al  pha ) 

equivalence  ( wmact,  Ymact),  (rrmnom,  vmnoin' 

equivalence  (  f f •  f).  (  dd,  d),  (xx,  x),  (anull*  if) 

call  second( timel) 

set  the  r  ow  dimensions 


1 da=250 
1  d  .i  f =  1 50 
1  dll  1=30 
1 dv=30 
Ida 1=25 
1 da0=30 
1  dt<0=30 

set  precision  of  solution  i e  01  =  IX 
prec  =  0000001 


Input  section 

read  in  from  standard  input 


58 


55 

60 

65 

70 

75 

80 

85 

90 

95 

100 

c 

98 


c 

o 

c 

101 

102 

St 

103 

107 

110 

111 

112 

c 


o 

c 

c 


read<5>  *)  n,  -9.  p»  i  n>  ou 
read<5>*)  <rnom<  i  ) »  i  =  l;  n> 
read<5/*)  <ract<i).i  =  l>n) 
do  55  i=l/ n 

r  04  d  <  5i  * )  <  1 1 1  (  i ,  .i ) ,  i=\,  n ) 
do  60  ,i  =  l ,  i  n 

read  <5.  *)  <  112<  i,  J  ) ,  i  =  l ,  n) 
do  65  i=l . ou 

read<5>*)<121(i.  ,i ) .  j  =  l#  n) 
do  70  j=l . i n 

read<5»  *)  <  122<  i»  j ) »  i  =  l<  ou) 
do  75  j  =  1 »  o u 

read(5>*H121r(i>  i ) .  i  =  1 ,  n ) 
do  80  j=l»  p 

r  ea d <  5>  * )  <  v <  i ,  ) ,  i  =  1 »  n ) 

read <5.  *) <  zpow< i ) »  i  =  l,  n) 
read<5*  *)  <  ss<  i  )  ^  i  =  l,  2*q) 
do  85  i=l . i n 

r ead<5>  *) <  uu<  2*i-l ,  J  ) ,  uu<  2*i,  .,i ) ,  i  =  l ,  q) 
do  90  i*l>  ou 

r  ead<  5;  * )  <  wmnom  ( 2*i-l ,  j) ,  YYmnom<  2*i <  j  ) ,  -i  =  1 .  q) 
do  95  i  =  l> p 

r ead<  5, *) <  aal pha<  2*i-l >  j  ) »  aal  pha<  2*i  >  i  )  >  .i  =  l .  q) 
do  100  i  =  li  ou 

r  ea  d  <  5,  * )  <  YYma  c  t  <  2* i  - 1  •  • ) »  YYma  c  t,  < 2*  i  >  j  ) .  ,i  =  1 ,  q  > 

i tmax=200 
do  98  1  =  1.  p 
read<5>  *»  end=101 ) 
naf=l 
error=.  15 

read<5,  *»  end=101)  naf»  error 


print  a  program  heading 
wr  i te<  6/  102) 

f ormat( //20x>  'Faul t  Diagnosis  Prosram— Tabl eau  Method” 
25x» ' (I nter polate  to  find  lambda)”) 
wr i te<  6»  103) 

format</25x,  " Nominal  Parameter s' ) 

wr i t  e  <  6 • 107)  (  'r  <  ' , i ,  r  n  om  < i ) »  i  =  1 »  n ) 

f ormat<4(3x.  a2>  i2,  a2»  elO.  4 ) ) 

Mjri  te(6.  110) 

f  or ma t < /25x >  'Test  Frequencies') 
do  111  i*l,  q 

wri  te<  6*  112)  i ,  s<  i  ) .  <  uu<  -i  •  i  ) ,  ,i  =  l ,  2*i  n) 
formats  1  Ox ,  '  s(  ",  i2,  .  ell.  4,  •'  +  ,i ' .  ell  4, 

2x,  'u='.  10<  f4  1,  '  +.('•,  f 4.  1,  5x)  ) 

wri  te<6>  ■*> '  #  of  allotnable  faults.  '  >  naf 

wr  i te<  6»  * ) '  all  owa  b  1  e  err  or :  x  .  err  or 

Compute  the  vector  bO  =  121r <vmact-122*u> 

do  130  i-1, q 
do  125  j=l . ou 
tmp<  j  )  ®cmplx<0  *0  > 
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120 

125 


130 

c 

0 

c 

c 


201 

204 

c 

c 

CXXXXX 


280 


290 


c 

295 

c 


S< 

300 


c 

c 

c 

305 


c 


fy 


C 


do  120  k=l . in 

t,mp<  j)=tmp<  i)+122(.i,  k)*u(k.  i) 
tmpC  j  )=Ymact<  .i,  i  )-tmp(  ..i ) 
do  130  J=l# n 
bO<i>  i)=cmplx<0.  »  0.  ) 
do  130  k=l» ou 

bO<  j »  i  )  =  bO<  i>  i  )+121r<  i»  k)#tmp<  k) 


Compute  the  vector  aO 

do  204  i*l» g 
do  204  i  =  1 » n 
aO<J.  i)=omplx<0.  .  0.  ) 
do  201  k=l .  n 

aO<  i,  i  )=aO<  i,  i  )  +  lll  (  ,i .  k)*bO<  k,  i  ) 
do  204  k=l.  in 

aO  <  i .  i  )  =aO  <  J  •  i  )  + 1 1 2  <  J  •  k )  * u (  k»  i  ) 


:  x.  x  x  x  x  x  x  x  x  x  x  x  x  x:  x  x  x  x  x.  x  x  x:  x  x  x  x  x  x  x  x  x:  x  x  x  x  x  x  x.  x  x  x  x  x  x  x  x  x  x  x  x  x  x  x  x  x  x.  x  x  x  x 

do  280  i=l ,  n 
rkthnm<  i )  =  1. 

dev< i )=<ract< i )-rnom< i ) >*100  /rnom( i ) 
i ter=0 
i t  n  ex  t=0 
xlam=-l . 

wr i te< 6. ' < /// ) " ) 
i ter=i ter+1 

i  f  <  i  t  er.  g  t.  i  tmax )  the  n 

write<6.*>"  iteration  limit  exceeded" 

go  to  900 

endi  f 

do  295  i  =  l,  r. 

r  k  t  h  <  i )  =r  k  t  h  nm  <  i )  *r  n  om  <  i  ) 
eva  1  ua  t  e  f 
do  300  i  =  l,  3 

cal  1  guadf  <  f  <  i*n-n+l )  •  rkth,  al  pha(  1 .  i  ) .  s<  i ) .  2pow»  111.  1  dl  1 1 . 
m.  1  dv.  aO  ( 1 .  i )  #  bO  (1.  i  ) .  n.  p.  tm  F’ ) 

conti nue 

f norm- norms g(  ff.  2*n*g) 

mr i t  e ( 6  >  * )  "  I TERAT I ON  1 •  i t  er .  "  I F I " .  f n  or m 

wr  i  te<  6.  * ) 

if<fnorm.  le.  1  e-10)  so  to  500 


stop  because  the  algorithm  appears  stuck  in  a  relative  min 
if(xlam.  eg.  0  )  wri  te<6<  305) 

format ( lOx,  37hcausht  in  a  point  of  relative  minimum) 
if<xlam,  es.  0.  >  stop 

call  iacobtif.  1  d-i  f  •  r  kth.  al  pha.  1  da  1 .  s.  n.  p.  g.  zpow.  111.  ldlll- 
v>  Idv.  aO.  IdaO.  bO.  ldbO.  tmp.  aa) 


c 


call  cxtorl<a.  Ida.  if.  ldif.  n.  f.  g) 
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c  par 

c 


315 
c  Far 

o 

o 

cddd 


cddd 

c 

c 


31? 

321 

c 

c 


& 

333 

32? 

331 

c 

33? 

c 

c 

c 

410 


420 


c 


455 

456 

457 


since  Jacob  computes  C  partial  f  /  partial  rfcth  3  the  last  n  column 
the  matrix  a  must  be  adjusted  to  set  i  partial  f  /  partial  rkthnm  3 
do  315  i  =  l,2*n*s 
do  315  J  =  1 ,  n 

a  <  i .  2*  p*  J )  =a  ( i  ,  2*  p*  s+  J  )  *r  n om  (  * ) 

use  libs  f  t  o  find  the  s  ear-  c  h  dire c  t  i  o  n 
llbsf  is  from  the  imsl  library 

n  f=l 
i  nd=l 
c  <  1 )  =0. 
c< 2 )  =  1.  e-5 
c  <  3 )  =  1 . 

call  llbsf<a.  1  da,  2*n*s.  2*p*s+n.  f.  Ida.  nf,  ind,  c»  d.  Ida.  ip.  ..if.  ieri 

NOTE:  the  r  part  of  d  is  normalised 

do  31?  J  =  l,  s 
do  31?  i-1 ,  F 

x  <  (  J  - 1 )  *  f  +  i  )  =a  1  f-  ha  <  i »  J  )  -0.  0*  d  <  <  J  - 1 )  *  f-+  i  ) 
do  321  i-1,  n 

xx  <  2*  f*  s+ i )  =  ( r k t  h nm ( i ) -0.  0*  d d ( 2* p*  s+ i ) ) *r n om ( i ) 
the  -0.0  is  the  lower  end  of  the  interpolation  data 
compute  function  If  1**2  at  lambda-0  -  25  -  5  -  75  -1 
do  33?  i lam=l . 5 
do  333  i-1,  s 

ca  11  s ua d f <  f ( i * n— n+ 1 ) ,  xx ( 2* p* s+ 1 ) , x ( ( i  —  1 ) * p+ 1 ) .  s(i).  z p ow • 

111.  ldlll.  m.  1  dv,  aO(  1 .  i  ) .  bO<  1 .  i  ) .  n.  p,  tmp) 

continue 

3<  i  lam)  =normss(  ff,  2*n*s) 
do  32?  J-l,  s 
do  32?  i-1 . f 

x< (  J-l '*p+i ) =x< ( J  — 1 )*p+i ) -  25*d( ( J-l )*F  +  i ) 
do  331  i-1.  r. 

xx  <  2*p*s+i  )  -xx (  2*f*s+  i  )  -  25#  dd<  2*p*s+i  )  *r  nom  (  i  ) 

the  -  25  is  the  increment  between  the  interpolation  data 

continue 

xlam-lambdat  s) 

update  the  alpha's  and  r 

do  410  J-l,  s 
do  410  i-1. p 

a 1 p  ha ( i .  J ) =a 1 p  ha  < i ,  J  >  +x  1 am*  d  < ( J - 1 ) *  p+  i  ) 
convr s=  true, 
do  420  i-1,  n 

i  f  <abs<  xlam*dd(  2*p*s+i  l  ) .  st.  prec)  Conors—.  false. 
rkthnm< i ) -rkthnm* i ) +xlam*dd< 2*F*s+i ) 
if < convr s)  so  to  500 

write  the  parameters  as  a  prosress  report 
write*  6.*)  rank  of  a  :  •' .  i  nt(  c<  4)  ) 

wr i  t  e  <  6 . 455 )  i t  er , x 1 am 

f  orma  t  < 1 Ox , 26  h n  or ma 1 i 2  e  d  r  at  it  era t i 0  n  .  i 3 , 2x . 

3h<  lambda-.  el2  5,  1  h)  ) 
wr i t  e ( 6 • 457 )  ( r k  t h  nm < i > . i  =  1 .  n ) 
format(5x, 4el5  5) 


write (6,  '<///)  •' ) 
so  to  290 


c  at  this  point  the  iteration  has  found  a  "feasible  solution" 

500  if<itnext.  e a.  iter)  so  to  900 

do  510  i  =  li  n 

510  r  k  t  h  (  i  )  =r  k  t  h  nm  <  i  >  *r  n  om  (  i  ) 

call  .iacob< -if.  1  d j f ,  rkth,  al F-ha,  ldal,  s,  n.  p,  s,  2pow<  111,  ldlll, 

&  Vi  ldv,  aO,  ldaOi  bO»  ldbO,  tmp,  aa) 

call  cxtorKa,  Ida,  jf,  1  d  >  f ,  n,  p,  s) 
do  515  i=l , 2*n*s 
do  515  j*l,  n 


515  a  < i ,  2*  p*  s+ i )  =a  (  i ,  2*  p#  s+  .i )  *r  n  om  (  i ) 

call  nulla(a,  1  da,  2*n*s,  2*p*<i+n,  anul  1 ,  Ida.  nldim,  dd,  tmp) 
c  t  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  t 

call  f  i  nder  <  anul  1  (  2*p*<«+1  ,  1).  Ida,  n,  nl  dim,  r  kthnm.  naf.  error-  beta) 
c  compute  anull*beta  and  put  it  in  d 

do  600  i=l, 2*p*s+n 


dd<  i )  =0. 

do  600  j*l,  nldim 

£.00  d  d  <  i )  =  d  d  (  i  )  +a  n  u  1 1  <  i ,  •• )  *  b  e  ta  (  • ) 


d  o  6 10  i  =  1  > 
do  610  i  =  l,  p 

£>10  a  1  p  ha  (  i ,  ..i )  =a  1  p  ha  <  i ,  .i )  +  d  (  <  J  - 1 )  *  f  +  i  ) 

convrs=.  true, 
d  o  620  i  =  1 ,  n 

i f < abs< dd< 2*p#s+i ) > .  at.  prec)  convrs=,  fal  se. 

£.20  rkthnm<  i  )=rkthnm<  i  > +  dd<  2*p<M+i ) 
i tnext=i ter+1 

wri te< 6, 630)  iter.  beta( nl dim+1 ) 
o-  beta  (  nl  dim+1 )  contains  the  best  norm  for  r 

£.30  f  orma  t  ( 1  Ox ,  n  orma  1  i  2  e  d  r  at  it  era  t  i  0  n  "  •  i  3,  2x  • 

S <  '(null  space  step:  rnorm=’.  elO  4.')  ') 

if(.  not.  convrs)  so  to  456 


0 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1  e  s  1 1 


0  wr i t  e  0 u  t  final  va lues 

?00  call  second( time2) 

wr i t  e ( 6 • 923 )  i t er- 1 ,  t i m  e2- 1 i m  e 1 

?23  forma t< /5x,  i5,  2x,  16h  iterations  and  ,  el 2.  6,  14h  sec  required) 

write (6,  925) 

?25  f orma t < /5x , 36 hParam e t er  Values  at  S o 1 u t i 0 n  P 0 i n t 

&  1 9 h%  dev.  from  nominal,  17h  error  in  sol  > 

do  950  i=l , n 

927  formats  lOx.  2hr<  .  i2,  4h)  =  .  ell.  4,  18x,  flO  2,  5x,  flO  2> 

?50  wr i t  e  < 6 , 927 )  i • r  k t  h  ran  ( i ) *r  n  om  < i ) ,  d  ev ( i ) , 

St  1 00  *  ( ra  c  t  (  i  )  -r  k  t  hnm  (  i  )  *r  n  om  (  i  )  )  /ra  c  t  (  i ) 

s  t  o  P 


end 


s u  br  0  u  t i n  e  find  er <  a  nu 1 1 ,  1  da  •  n .  nldim,  r  k  t h nm ,  na  f •  err .  b  e  ta ) 
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integer  negl,  nldim,  n»  1  dr,  1  eg  1  ( 30 ) .  inegi<30).  imK30l . 

St  i  S4M  e  <  30 ) .  i  n  d  ex  <  30 ) 

real  rkthnm(l).  diff<30),  rdiff<30),  rnul 1 < 30. 30) .  r<30), 
h.  bsaue(30).  anul  1  (  1  da,  1 ) .  beta(l).  iijk<2500)*  c<  4) 

c 

1  dr  =30 
err or =  err 
do  30  i  =  l  *  r. 

30  d  i  f  f  (  i  )  =  1 .  -r  k  t,  h  nm  (  i  ) 

c 

c  determine  which  parameters  are  already  at  nominal 

i  i=0 

do  1 00  i = 1 .  n 

i  f  (abs<  di  f  f  <  i  )).  le  error)  then 
i i=i i+1 
i  egl  < i i >  =  i 

e  n  d  i  f 

1 00  continue 

c 

n  e  g  1  =  i  i 
nonegl=n-negl 
c  setinegl 

130  ii=0 

do  150  i  =  l,  n 
do  140  ,i  =  l  .  negl 

140  i  f  <  i.  eg.  i  egl  ( -• ) )  so  to  150 

i i=i i+1 
inegl( ii)=i 
150  continue 

wr i te< 6, *)  '  possible  faults"*  < i negl < i ) ■  i  =  l  •  i  i  ) 

c 

sumold=l  e+10 
cnnn  icnold=n 

c 

i  ..i  ob=0 

c  add  a  sroup  of  the  possibly  bad  to  the  good  < nadd=nonegl-naf ) 

na  d  d=  n  o  n  e  g  1  -  na  f 
if (nadd  It  0>  nadd=0 

160  call  incr<  index*  nadd*  no  negl,  i  job) 

0322  wri  te(  6,  *)  '  i  ndex-' .  (index(i),  i  =  l»  nadd) 

if<i  .iob  eg.  0)  so  to  240 
c 

do  170  i=l,  nadd 

170  i e g 1 <  n e g 1  +  i )  =  i ne g 1 < i ndex ( i ) ) 

n  e  g  n  s=  n  e  g  1  +  na  d  d 
do  ISO  i=l .  negns 
rdiff(i)=diff(iegl(i)) 
do  180  ..i  =  1 ,  nidi m 

180  r  n  u  1 1  (  i .  ,i )  =a  n  u  1 1  (  i  e  g  1  <  i  )  •  i ) 

c 

i  nd=  1 
ncol s=l 
c  < 1 ) =0 
c  ( 2  >  =  1 .  e-5 
c  <  3 )  =  1 

call  1 1  t*gf  <  rnul  1 .  I  dr,  negns,  nl  dim,  r  di  f  f  *  ldr*  ncols*  ind*  c* 
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beta.  1  dr,  iwk,  wk,  k) 


compute  a  trial  value  of  r 
do  190  i  =  l,  n 
r  < i )=rkthnm< i ) 
do  190  j  =  l,  nldim 
r < i )=r< i )+anul 1 ( i ,  i  )*beta<  J ) 


c  compute  a  performance  index  to  evaluate  the  fault  combination 

s  um=0. 

do  210  i  =  l,  nems 
err xx=a  b  s  <  r  < i e 3 1 < i ) ) - 1 .  ) 
s  um= s  um+  er r xx**2 
210  continue 

0  2  2  2  2  wr  i  t  e  <  6 ,  * )  •'  s  um :  ' ,  s  um 

sum=sum**.  5 

cxxxx  wr  i  t  e  <  6,  * )  y  i  e  e  1  ' .  <  i  e  3  1  (  i  ) ,  i  =  1 ,  n e 3 n s ) ,  '  s um :  " ,  s um 

c 

OMV  M  M  K>  M  M  U  MM  M  *-.*  MM  M  MM 

c  r  e  e  c  t  a  s  o  1  u  t  i  o  n  with  n  e  sa  t  i  v  e  par  am  e  t  er  s 

do  215  i=l,  n 

if<r<  i).  It.  0.  )  so  to  160 
215  continue 

cvvvvvvvvvvvvvvvvv 
c  compare  to  all  previous 

if  (sum.  It.  sumold)  then 
s  um  o  1  d=  s  um 
do  220  i=l, ne^ns 
220  i save< i l=i e^l < i ) 

do  230  i=l, nl dim 
230  bsave( i )=beta< i ) 

endi  f 
3 o  to  1 60 
240  continue 

d  o  250  i  =  1 .  nldim 

250  b  e  ta  < i )  =  b sa v e  < i ) 

beta<  nl dim+1 ) =sumol d 

r  e  t  ur  n 

end 


subroutine  inert  index,  ni>  n,  i  iob) 

This  s  u  br  o  u  t  i  n  e  in  or  em  e  n  t  s  the  i  n  d  ex  v  e  c  t  or .  i  n  d  ex. 

The  vector  index  should  list  the  indices  of  n  items  taken  ni  at 
a  time. 


variable  list 


i  n  d  ex  i  n  p  u  t/  o  u  t  put:  The  i  n  d  ex  v  e  c  t  or 


i  -i  ob 


i  nput 


i  n  p  u  t 


Order  of  index 


Or d  er  o  f  full  set 


i nput/ output:  Control  on  input  0=initiali2e 

1 = i ncrement 
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Flag  on  output 


0=  e  n  d 


integer  index(ni).  ni,  n>  iiob 

if<(ni.  eg.  0)  and.  (iJob.  e-a  0))  then 
i j  ob=l 
r eturn 

else  if((ni.  e<i.  0).  and.  (Uob.  ei.  D)  then 
i  .j  ob=0 
r  e  t  ur  n 

else  i  f  <  i  j  ob.  eg.  0)  then 
do  110  i  =  l ,  ni 
i ndex( i )=i 
i i ob=l 
r  e  t  ur  n 

else  i  f  <  n  i .  ei.  n )  then 
i j  ob=0 
return 


el  se 


endi  f 
end 


i ndex< ni )=i ndex< ni ) +1 
do  130  i  =  ni . 2, -1 
if(index(i).  at.  n-ni  +  i)  then 

i n  d  ex ( i  - 1 )  =  i  n  d  ex  < i-l)  +  l 
i ndex< i )  =  i ndex( i-1 ) +1 

endi  f 
conti nue 
do  150  i=2>  ni 

if(index<i).  st.  n-ni  +  i  >  index!  i  )  =  index(  i-1  )  +  l 

if  <  index!  1).  gt.  n-ni  +  1)  i,iob=0 

return 


subroutine  nulla! a,  Ida,  nrows,  ncol  s>  anul  1  >  ldnull.  nldim,  s, 

This  s u br o u t i ne  c om p u t e s  the  nulls pa c e  o f  the  ma tr i x  a 
Initially  the  null  space  dimension  is  set  to  ncol s-nr ows 
if  this  is  positive  or  2  er  o  o  t  h  erw  i  s  e.  The  r  o  u  tine  then 
checks  the  singular  values.  The  null  space  dimension  is 
incremented  for  ev  erv  singular  va  1  u  e  hi  h  i  c  h  is  too  small 
relative  to  the  larsest< first) 


var iable  list 

a  input:  matrix  for  which  null  space  is  desired 

a  is  d e s tr ov e d 


input  row  dimension  of  a 


nr  ow  s 


input  number  of  rows  in  a 


ncol  s 


input:  number  of  columns  in  a 


i 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


10 

20 


anull  output:  the  columns  of  v  span  null  Cal 
ldnull  input:  row  dimension  of  m 

nldim  output:  dimension  of  nullCal 

s  work:  work  array  of  dimension  mi  n<  nrows#  ncol  s) 

containing  the  singular  values 

wk  work:  work  array  of  dimension  2*max( nrows* ncol 


integer  Ida,  ldnull,  nrows,  ncols,  nldim 
real  a<  1  da,  1 ) ,  anul  1  <  1  dnul  1,  1 ) »  s<  1 ) ,  wk(  1 ) 

call  lsvdf  <a*  Ida,  nrows,  ncols,  b,  0,  0,  s,  w.k,  ier) 

nidi m= ncols-  nr  ow  s 

if<nldim.  It.  0)  nldim=0 

do  10  i  =  l ,  mi  n0<  nrows,  ncol  s) 

if<s<i)  le.  1.  e-5*s<l))  nl  dim=nl  dim+1 

do  20  i=l, ncol s 

do  20  4  =  1,  nldim 

anul  1  <  i  >  ..i )  =a <  i ,  ncol  s-nl  dim+  i ) 

return 

end 


The  following  subroutine  listed  in  Appendix  A  are  also  required 

guadf 
..i  a  c  o  b 
cx  t  or  1 
n  or  m  s  ■=» 
lambda 
second 


